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PREFACE 


An Improved Error Assessment for the GEM-T1 
Gravitational Model 

Several tests have been designed to determine the correct error 
variances for the GEM-T1 gravitational solution which was derived 
exclusively from satellite tracking data. The basic method employs both 
wholly independent and dependent subset data solutions and produces a 
full field coefficient by coefficient estimate of the model uncertain- 
ties. The GEM-T1 errors have been further analyzed using a method based 
upon eigenvalue-eigenvector analysis which calibrates the entire 
covariance matrix. Dependent satellite and independent altimetric and 
surface gravity data sets, as well as independent satellite deep 
resonance information, confirm essentially the same error assessment. 
These calibrations (utilizing each of the major data subsets within the 
solution) yield very stable calibration factors which vary by 
approximately 10$ over the range of tests employed. Measurements of 
gravity anomalies obtained from altimetry were also used directly as 
observations to show that GEM-T1 is calibrated. Based upon these 
calibrated error estimates, GEM-T1 is a significantly improved solution 
which to degree and order 8 is twice as accurate as earlier satellite 
derived models. By being complete to degree and order 36, GEM-T1 is 
much larger than earlier gravitational solutions calculated from direct 
satellite tracking and has significantly reduced aliasing effects that 
were present in previous models. The mathematical representation of the 
covariance error in the presence of unmodeled systematic error effects 
in the data is analyzed and an optimum weighting technique is developed 
for these conditions. This technique yields an internal self- 
calibration of the error model, a process which GEM-T1 is shown to 
approximate. This geopotential field with calibrated error estimates, 
predicts 25 cm for the radial RMS uncertainty of the TOPEX orbit. The 
TOPEX Mission has a requirement for 10 cm radial orbital modeling which 
is needed to support the oceanographic applications of a high quality 
spaceborne altimeter. 
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I . INTRODUCTION 


The least-squares method yields a formal estimate of a solution's 
accuracy. This estimate is generally optimistic due to imperfections in 
the mathematical model employed and suspected, yet unmodeled, systematic 
error sources resulting in non-random observation residuals. Hence some 
form of verification and calibration of errors is required to reliably 
establish the "true" accuracy of the results. This verification is made 
difficult in circumstances where the solution output is a mathematical 
model of a physical phenomenon whose empirical coefficients taken indi- 
vidually are not directly observed. This situation applies here to the 
problem of estimating the accuracy of a geopotential field determined 
from satellite tracking data. Such models have grown in size and 
resolution as tracking instrumentation has advanced and many satellite 
orbits became available. It is desirable to validate such a solution 
completely and even though recent models now contain more than 1300 
individual coefficients, this has become feasible with the main method 
employed here (Lerch, 1985) and the .use of high speed vectorized 
computers. It is also desirable and feasible to validate the entire 
error covariance matrix through an eigenvector analysis as developed in 
this report. 

We discuss here the error calibration of "satellite only" Goddard 
Earth Models (GEM). These fields have been computed from restricted 
observation sources sensitive only to satellite dynamics. To serve the 
objectives of error calibration, data rich in high frequency local 
information have been excluded when deriving these fields. An accuracy 
assessment of these models benefits from a direct comparison with such 
information derived from satellite altimetry and surface gravimetry 
which both independently describe the gravitational field over vast 
regions of the earth's surface with a resolution well beyond the limits 
of satellite dynamical sensitivity. Tests of the geopotential solution 
with global sets of these independent local observations have yielded 
important statistics on the true accuracy of the gravitational models. 
Historically, the Goddard Earth Models have been calibrated through the 
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use of these and other independent data as described in Lerch (1985), 
Lerch et al. (1974, 1979, 1985, and 1986), Wagner (1983), and Wagner and 
Lerch (1978). These calibrations were made using a variety of tests and 
comparisons including those already mentioned with surface gravimetry 
and altimetry, and also utilizing deeply resonant orbits and assessments 
of orbit solutions when fit to tracking data excluded from the 
solution. All such tests in the past used data sets that were 
essentially independent of the model. 

For example, an earlier validation method, and one which has had 
some success, was to divide the final solution into independent pieces. 
Independent solutions were computed from these highly decimated portions 
of the original data set to assess solution weighting factors and scaled 
formal statistics. In the past, we have used this approach to compare 
the laser vs. optical contributions in the GEM models (Lerch et al . , 
1985). But this type of test has its limitations, for these subset 
solutions seldom approach the accuracy of the complete original field. 
Nonetheless, calibrations using this method have proven beneficial in 
establishing rough estimates of field accuracy while also testing data 
weighting approaches and locating suspect data contributions. 

The fundamental objective of our analysis in the GEM models is to 
develop the most accurate and well resolved spherical harmonic expansion 
of the gravitational field, and at the same time, obtain reliable esti- 
mates of its uncertainty. As the work progresses, the information 
included in the models is upgraded and made more complete through the 
incorporation of additional data types. Independent data tests are more 
difficult to construct as these formerly independent observations are 
merged into the solution. Obviously, the most accurate and complete 
models require the best observations to resolve the geopotential. The 
final solution inevitably contains all of the strongest data sets. 
Therefore, independent data are increasingly more difficult to find and 
the desired accuracy of the model outweighs the needs of specific error 
calibrations which require the exclusion of independent test data sub- 
sets. Our latest GEM-T1 solution (Marsh et al., 1988) has been developed 


2 


to serve as the base model for a whole new series of comprehensive 
geopotential solutions. Future solutions will contain virtually all 
significant satellite tracking data sets, and combination models will 
also contain surface gravimetry and satellite altimetry. 

In anticipation of the need to verify the accuracy of increasingly 
more comprehensive and precise gravitational models, a new technique was 
sought to achieve reliable full field calibrations. This new method 
developed by Lerch (1985) uses major dependent subsets of the full 
solution and is described further in the report on GEM-T1 (Marsh et al . , 
1987). This method of calibrating a solution is the main focus of this 
paper. We will show that in addition to the verification of the model 
uncertainties, this technique yields valuable information on the proper 
weighting of the data subsets comprising the field, and on the proper 
relative weighting of new data when they are introduced into existing 
models . 

This calibration technique is based upon the comparison of the 
complete model with a nearly equivalent one lacking the specific contri- 
bution of an individual observation subset. In this case, the subset 
field is not substantially reduced in performance and the complete 
complementary data set and the original model undergo simultaneous 
calibration. Since an individual data subset of GEM-T1 does not signif- 
icantly span all of the eigenvectors of GEM-T1 , it is important to test 
the calibration on each of the major data subsets individually. Hence, 
this report shows the error assessment of GEM-T1 using all of the major 
subsets comprising the solution. Furthermore, since GEM-T1 lacks gravi- 
metry, altimetry and deep resonant satellite information, the calibra- 
tions obtained through this analysis are compared with that obtained 
through more traditional calibrations using independent data. Finally, 
given that the solution uncertainties (error variances) alone of a 
satellite-derived model are insufficient to predict the performance of 
the field on specific orbits, a more complete calibration and validation 
of the error covariance matrix of the geopotential solution was 
required. This analysis was extended to the calibration of the 
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eigenvectors of the GEM-T1 error covariance matrix which is a major 
element developed within this report. 

Since the primary method of calibration subtracts a pair of 
solutions with and without a given subset of data, there is concern that 
a bias due to some unmodeled effect may exist in each of the solutions 
which would be eliminated through the subtraction and therefore remain 
undetected. In order to account for this type of problem, calibrations 
were also made which test the gravitational solutions directly with 
independent surface based and deep resonant observations. Gravity 
anomalies obtained from altimeter data and deep satellite resonant 
passage observations were employed for this purpose. 

A principal interest in these error assessments is to predict the 
performance of the geopotential on the yet-to-be-flown TOPEX satel- 
lite. Gravitational modeling error is the dominant error source which 
can degrade the precise ephemerides sought for this mission. An 
accuracy level of 10cm in the radial component of the TOPEX orbit is 
required which is substantially beyond the capabilities of any existing 
gravitational models. GEM-T1 is a preliminary model developed for the 
TOPEX Project and has reduced the field uncertainty for this mission. 
An assessment of TOPEX orbit modeling capabilities will also be pre- 
sented based upon the calibration of the full GEM-T1 covariance matrix 
since the covariance terms significantly affect this error estimate. 

An Appendix in three sections (A1, A2, A3) is presented to 
elaborate on the estimation theory employed in GEM-T1 and in this 
report. Because of the systematic error effects in the data and other 
unmodeled force effects, the mathematical validity of the representation 
of the covariance error is analyzed in Appendix A. Also included in the 
Appendix is an optimum process for simultaneously converging on all data 
subset weights. These weights compensate for the unmodeled systematics 
and provide an automatic self-calibration of the error covariance 
internal to the solution. It is shown that the method of solution for 
GEM-T1 approximates this process. 
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II. GEM-T1 CALIBRATION OVERVIEW 


GEM-T1 was derived from about 800,000 tracking observations on 17 
satellites (Table 1) consisting of laser, optical and Doppler data. 
When combining data from different instruments and different orbits 
whose sensitivity to non-gravitational signals vary widely, strategies 
must be developed and tested to obtain optimal data weighting factors. 
And to account for the general problem of optimistic formal statistics, 
we have found it necessary to scale down considerably the entire data 
set of the solution in order to yield a properly calibrated solution. 

The solution is obtained by simultaneously solving for a 36x36 
gravitational model, 66 dynamic tidal terms having long period orbital 
perturbations and earth orientation parameters for the post- 1979 time 
period. A total of nearly 2500 unknowns comprised the solution after 
back-substitution of approximately 7000 satellite-arc-dependent 
parameters (e.g., orbit elements, drag and solar radiation pressure 
coefficients and doppler biases). 

The post-fit performance of a solution on the tracking data — the 
so called "residual analyses" — plays a significant role in the design 
of the solution, the selection of data weights, and in the final error 
assessments. Marsh et al., (1987) discuss this problem and show in 
considerable detail relatively large orbit residuals for laser systems 
as compared to their precision level. In Appendix A2 these tracking 
residuals are analyzed on a pass- by- pass basis for some of the laser 
systems and show considerable biases which are unmodeled in GEM-T1. 
Hence, this verifies the need for the downweighting factors to 
compensate for unmodeled biases and to test the calibration of the error 
covariances. 

This report utilizes all of the major data sets and will more 
completely verify the GEM-T1 accuracy assessments through a 
comprehensive calibration of the solution variances and through a new 
technique utilizing eigenvector analyses which will test the entire 
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TABLE 1 


SATELLITE DATA 


SEMI MAJOR 


SATELLITE 

1 LAGEOS 

2 STARLETTE 

3 GEOS-3 

4 PEOLE 

5 BE-C 

6 GEOS- 1 

7 GEOS-2 

8 Dl-C 

I 

9 DI-D 
TO SEAS AT 

1 1 OSCAR-T4 

12 ANNA-1B 

13 BE-B 

14 COURIER-1B 

15 TELSTAR-1 

16 VANGUARD-2RB 

17 VANGUARD-2 


AXIS {km.) 

ECC 

DEG 

TYPE 

12273. 

.0038 

109.85 

LASER 

7331. 

.0204 

49.80 

LASER 

7226. 

.0008 

114.98 

LASER 

7006. 

.0164 

15.01 

LASER 

7507. 

.0257 

41.19 

LASER 




CAMERA 

8075. 

.0719 

59.39 

LASER 




CAMERA 

7711. 

.0330 

105.79 

LASER 




CAMERA 

7341. 

.0532 

39.97 

LASER 




CAMERA 

7622. 

.0848 

39.46 

LASER 




CAMERA 

7170. 

.0021 

108.02 

LASER 




DOPPLER 

7440. 

.0029 

89.27 

DOPPLER 

7501. 

.0082 

50.12 

CAMERA 

7354. 

.0135 

79.69 

CAMERA 

7469. 

.0161 

28.31 

CAMERA 

9669. 

.2429 

44.79 

CAMERA 

8496. 

. 1832 

32.92 

CAMERA 

8298. 

.1641 

32.89 

CAMERA 


SIGMA {a)={±) 


IN GEM-T1 


# OF 
ARCS 

57 

46 

36 

6 

39 

50 

48 

43 

28 

46 


14 

14 

13 

30 

20 

10 

30 

10 

10 


# OF 
OBS 

RMS 

RESID. 

S 1 GMA* 
WEIGHTS 

144527 

10cm. 

112cm. 

57356 

20cm. 

224cm. 

42407 

70cm. 

816cm. 

4113 

90cm. 

816cm. 

64240 

50cm. 

577cm. 

7501 

2 arcsec 

5.6 arcsec 

71287 

70cm. 

667cm. 

60750 

1 arcsec 

8.9 arcsec 

26613 

80cm. 

816cm. 

61403 

1 arcsec 

8.9 arcsec 

7455 

150cm. 

816cm. 

2712 

2 arcsec 

7.3 arcsec 

11487 

100cm. 

816cm. 

6111 

2 arcsec 

8.9 arcsec 

14923 

70cm. 

707cm. 

138042 

.5cm/sec 

7cm/sec 

63098 

Icm/sec 

8cm/sec 

4463 

2 arcsec 

4.5 arcsec 

1739 

2 arcsec 

4.5 arcsec 

2476 

2 arcsec 

4.5 arcsec 

3962 

2 arcsec 

4.5 arcsec 

686 

2 arcsec 

4.5 arcsec 

1299 

2 arcsec 

4.5 arcsec 



covariance matrix. To do so convincingly, different accuracy estimates 
and error calibrations will be undertaken and compared spanning the 
spectrum of eigenvectors in GEM-T1 . The stability of results across 
this spectrum will give insight into realistic field accuracy and will 
demonstrate the reliability of this determination. The coefficient 
uncertainties for GEM-T1 are shown in Figure 1.1. They are also 
presented as an rms coefficient error by degree in Figure 1.2. 

A review of the solution technique (Section III) utilized to 
achieve the weighting factors employed in the computation of GEM-T1 is 
useful to the discussion of the calibration method which is presented in 
Section IV. 
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ESTIMATED ERRORS FOR GEM-T1 COEFFICIENTS 

Figure 1.1 


ERROR x 10 

b 

□ 

■ GEM-L2 
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terrestrial gravity field 


J 2 4 6 8 1 0 12 1 4 16 1 8 20 22 24 26 28 30 32 34 36 

DEGREE 

GEM-L2 — Previous gravity model (20 x 20 field) 

GEM-T1 — TOPEX gravity model (36 x 36 field) 

PGS-T2' — TOPEX model: GEM-T1 without low inclination satellites 

RMS Coefficient Error per Degree 
Figure 1 .2 
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III. TECHNIQUE AND DATA WEIGHTING FOR THE GEM-T1 SOLUTION 


The method of solution is a modified --least squares process which 
minimizes the sum (Q) of signal and noise as follows: 


Q = 


1 

fl, ,m 


l ,m 


+ S 


&,m 


r 

l 1 

t obs a 
i 


it 

2 

t 


where the signal is given by 


(3.1) 


C 2 . spherical harmonics comprising the solution 

S, f m > H,m* coefficients; and 



is rms of the coefficients of degree a (a priori 
rule) and is introduced to permit larger solu- 
tions to degree and order 36x36. This law, 
based upon Kaula's rule, has been obtained inde- 
pendently from studies of the spectra of the 
Earth's gravity field and is used here to repre- 
sent the observed power within the geopotential. 


and the noise by 


observation residual (observed-computed) 
for the ith observation of satellite 
tracking data set (type) t; and 


RMS of observation residuals (generally 
• significantly greater than a priori 
data precision) " ' 

downweighting factor to compensate for unmodeled 
error effects in the data (ideally f=l). 


The weighting factor f is used to scale the data contributions; 
unmodeled errors remain in the data and this scale factor is utilized to 
obtain realistic error estimates from the solution statistics (variance- 
covariance matrix). A proper combination of the data normals with the 
well-known size of the a priori signal (o ) is the result. Since f is 
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not well known, as it depends upon the biases in the data (see Appendix 
A1), it must be estimated in the development of the solution from the 
data itself. The value of f was obtained from an extensive calibration 
procedure which we describe in what follows. Presently we have found 
that the value f = .02 is required since residuals in the tracking data 
are systematic and highly correlated particularly within a given 
tracking pass. 

When minimizing Q in (3.1) using the least -squares method, the 
normal metric equation and error covariance is obtained as follows: 

Mz _ R are the normal equations, where z is the 

solution, R is the vector of residuals, and 

is the approximate form for the variance- 
V zz = N covariance error matrix which must be (3.2) 

calibrated. 

The complete expression for V zz is given in Appendix A2 which accounts 
for systematic unmodeled effects. It is shown therein that, with the 
proper downweighting factors, M ^ is a reasonable approximation. 

Certain characteristics of the solution are important. The 

2 

relative weighting factors (1/a.) are first obtained as the rms ( a, ) 

L t 

of the observation residuals from the solution for each satellite data 

type t. In converging o t , a few test arcs of each satellite data type 

are employed. A value of f is also iterated upon and is obtained 

principally by comparing the solution with independent gravity anomalies 

in the usual calibration procedure. For example, we have found (Marsh, 

et al., 1988) the overall weighting factor, f = .02, provides for 

realistic error estimates of the coefficients in GEM-T1 . It is shown 

also (Appendix A 1.2 and A 1 . 3 ) that this downweighting accounts for the 

systematic trends in the data residuals which are not random and are 

strongly correlated over a pass of tracking data. The residuals in 

post-fit analyses (Appendix A 1 . 1 ) are predominantly trended and 

typically fall on a straight line over time for a given tracking 

interval. If the residuals were random then f=1 (which is the ideal 

case) would be realized for the least -squares weighting. The need to 
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use f=.02 for the solution scaling factor is interpreted as an 
indication that non-gravitational signals must be better accommodated. 
Usually a few iterations of combining and fine-tuning the weights for 
the normals for each satellite data type are required to converge on 
proper weighting factors. In Appendix A3 an optimal technique is 
developed allowing f(f t ) to converge simultaneously on all subset data 
types. This was applied to a certain extent in the development of 
GEM-T1 as indicated by the weights given in Table 1. 

Another important aspect of the error estimation process is the 
effect of the a priori signal on the solution (in the minimization of Q 
in (3.1)). The signal consists of the gravitational coefficients with 
an a priori estimate of their size, cr , like the modified version of 
Kaula's rule (— 1 10 ) which is used here. The normal least 
squares adjustment solely minimizes the noise (observation residuals) to 

obtain the solution. The process used here is a modified least squares 
adjustment which minimizes both the signal (the size of the potential 
coefficients which are constrained with an a priori power spectrum) and 
the noise simultaneously. This latter approach (Lerch et al., 1977 and 
1985) permits a more satisfactory and accurate estimation of high degree 
and order terms. This process is also referred to as collocation 
(Moritz, 1978, 1980). The process is necessary for solution stability 
especially in the determination of the high degree portion of the 
field. While there is some sensitivity within the data to these high 
degree terms, there is a lack of separability for the entire set. 

With the normal least squares approach (noise-only minimization) 
there is a problem of separability due to the strong correlation between 
many of the high degree coefficients. The absence of collocation (GEM-T1 
without the Kaula constraint) results in excessively large power in the 
adjustment of the potential coefficients. Figure 2 illustrates the 
instability of the least-squares solution when collocation is not 
used. A satellite-derived gravity solution has been solved without 
collocation which is evaluated using a global set of independent gravity 
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Gravity Model Comparison with 1114 
5°x5° SEASAT Gravity Anomalies 

Figure 2 
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anomalies. An unrestricted high degree field performs poorly due to 
excessively large adjustment in the coefficients which is normally 
circumvented in the standard least-squares method by solving for a 
smaller sized field. Unfortunately, by restricting the size of the 
field, one also is requiring the higher degree terms above the field 
limits to be constrained absolutely to zero. Figure 2 also shows the 
disadvantage of this approach where the smaller sized field (PGS-3067) 
contains aliasing in its coefficients and does not perform well. (The 
abbreviation PGS stands for Preliminary Gravity Solution.) The aliasing 
signal sensed in the data above the field limits is absorbed into the 
adjustment of the lower degree coefficients. The best approach is seen 
with the least- squares collocation (or constrained) solution, GEM-T1, 
with a complete solution of a 36x36 field in harmonics. 

Since the a priori signal matrix contains only diagonal terms 
which are added to the data normals, it is interesting to compare the 
relative size of their contribution to the overall data normals. As 
seen in Table 2, the satellite normal equations have considerably larger 
diagonal terms than does the collocation matrix. At lowest degree, the 
ratio of the two contributions indicates that the collocation 
contribution is barely visible. It is only at highest degree that the 
collocation contributions become a large percent (30% ) of the diagonal 
terms. Since collocation stabilizes the entire solution, it does so by 
indirectly controlling ill-conditioning due to correlation within the 
system of observation equations. 
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Table 2 


Ratio of Diagonal Elements of 
GEM-T1 vs. Signal (Kaula) * Matrix 



ratio** of diagonal elements for specific terms 

Degree 

Order 

Sectorial 

(0 

m=1 

nW 

2 

400,000,000 

4,000,000,000 

5 

5,000,000 

1,700,000 

10 

400,000 

200,000 

15 

30,000 

30,000 

20 

1,500 

10 

25 

300 

7 

30 

50 

6 

35 

7 

3 


* 


Collocation/ Kaula Weight = 


1 

-5-1 2 

1 to 

,V2 i . 


- 2xlO'V 


** — 


DIAG 


GEM — T1 


DIAG 


Kaula J 
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IV. MATHEMATICAL DESCRIPTION OF ERROR CALIBRATION 
BASED UPON FIELD SOLUTIONS 


The most direct method of calibration entails the computation of a 
field solution wholly independent of the data used in GEM-T1 and a 
comparison of GEM-T1 with its results. The coefficient differences 
( AC am ) combined with the error estimator (e ) between GEM-T1 and this 
test solution will be used to compute the following quantity as a sample 
calibration factor for each spherical harmonic as follows: 


The ensemble of the sample will be utilized to determine a calibration 
factor for terms of degree a, or order m in GEM-T1 together with the test 
field from the formulae developed in this section. An overall 
calibration factor can then be obtained from the rms average of the 
sample calibration factors. 


GEM-T1 will also be compared to dependent solutions using the 
GEM-T1 data set with individual subsets eliminated. Calibrations of the 
same form as shown in (eq. 4.1) will be developed using these models. 
Since the estimated errors in both of these approaches will involve 
both GEM-T1 and the test field in a mutual calibration, it is important 
for reliability of the entire error spectrum that a significant number 
of test solutions based upon a variety of data sets be employed. 

The models which will be tested are complete to degree and order 
36 containing 1363 terms. A separate calibration factor for each 
harmonic is computed and for overall model assessments, these factors 
are statistically combined by both degree and order. 

A mathematical description of the error calibration is now 
presented with a derivation of its components. 
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First, let x and x be the respective errors in two gravitational 
models F and F consisting of coefficients C and C respectively and let 
E( ) define the expected value. It follows that: 


C 


= C 


true 


x 


C 


C, + 
true 


x 


and 


C 



= x - x 


represent the coefficient differences. Both x and x are errors within 
each solution of unknown magnitude. The covariance matrix of these 
errors is: 

V- — = E(x x T ) for solution F 

X X 


and similarly 

T 

V = E(x x ) is the covariance matrix for solution F. 

X X 

The cross covariance for x and x, non-zero when they share observations, 
is V x — . The covariance V for the difference in the two solutions is 
obtained as: 


V = E[(x-x) (x-x) T ] = E(x x T ~2x x T + x x T ) 


= V 

X X 


2V - 


x x 


+ 


V- - 

X X 


(4.2) 
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For solutions which are completely independent (i.e. they share no data 
subsets) as noted above 


V - 
x x 


0 


Therefore, from (4.2) we find 


V 


V 

x x 


+ V 


X X 


(4.3) 


For cases where the data in F is totally contained within the F 
solution, we will show that: 

V - = V 

XX XX 

Then, again using (4.2) the variance of the difference solution (C-C) 
is simply: 


The complete solution has (n + nl) observations while the subset 
solution only contains n observations. Let the contribution of the 
least - squares normal equations be denoted as Nx = R and Nix = R1 
respectively for the set of n and nl observations. Then for the 
combined set of observations n + nl: 

(N + N1)x = R + R1 (4.5) 

are the set of complete normal equations whose covariance for the case 
of random noise is: 

(N + Nl)' 1 = V xx (4.6) 
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Likewise, for the subset solution 

(N)x = R ( 4 , 7 ) 

with a covariance 

(N)_1 = V x x (4.8) 

See Appendix A2 where this is shown to be an approximation. 

From (4.5) and (4.7) 

xx T = N -1 R (R + R1) T (N + N1 )" 1 ( 4 . 9 ) 

and 

xx T = N _ 1 (RR T + RR1 T ) (N + N1 )" 1 (4.10) 

Since R and R1 are based upon independent observations they are 
uncorrelated 

E(R R1 T ) = 0 (4.11) 

The expected value for E(x x), using (4.10) and (4.11) is then 
E(x x T ) = V - = N " 1 N E( x x T ) N (N + KM ) -1 

X X 

= (N + N1 )" 1 = V xx (4.12) 


— — T _-j 

since E(x x ) = N . Thus we have the general result for V (the 
covariance of solution differences) 

y _ y + y for C and C being 

x x xx independent solutions, and 
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(4.13) 


V = V V 

x x 


X X 


for C with data wholly contained 
within the data of C 


Since these results (4.13) are approximations based upon random noise, 
they need to be calibrated which is done in the remainder of this 
report. To this end, we define several quantities used in our field 
calibrations : 

For two fields F and F: 

F : C n , S n , o's (coeff. errors) 

Am * Am’ 


F : C , S , o's 
Am 7 am 7 


RMS Jt (AF) = 


a AC„ + AS 

I , , 

m=0 2a + 1 


a,m 


1/2 


2 2 

i °(C ) + »(S. ) 

£ a,m a,m 

m=0 On 1 

2a + 1 


1/2 


(4.13a) 


since there are 2A+1 coefficients per degree. Further, for order m 


RMS (AF) 
m 


36 AC* + AS?, 
^ am am 

a=m 6 (36-m+1) 
m 


1/2 


6 = 1 , m-0 

m ’ 

= 2, m*0 


m 


2 2 
36 °(C ) + ^ ( S ft rr, ) 

v am am 

6 (36-m+1 ) 
a=m m 


1/2 
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since, for a 36x36 field, there are <5 m (36-m+1) coefficients per order. 
From the derivation given we know that the expected value of the average 
coefficient differences for a given degree, e^, is: 


el = E( RMS^F . 

When F is independent of F as in (eq. 4.2), the derivation (4.13) shows 
that 

2 2 -2 

e * = a i + (4.14) 

and when data in F is wholly contained within F as in (eq. 4.4), the 
derivation shows that 


e 


2 



(4.15) 


Based on the above equations and (eq. 4.1), the calibration factors we 
us® throughout the next section are defined by: 


k 


% 


rms a 



for degree fi, 


(4.16) 


k 


m 


rms 

m 


e 


m 


for order m 


(4.17) 


rms 


"a,m 


l ,m 


’£,m 


for an individual coefficient pair (4.18) 
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V. CALIBRATION OF GEM-T1 USING DEPENDENT SUBSET AND 
INDEPENDENT TEST SOLUTIONS 

5.0 Introduction 

Many investigations utilizing the gravity model also require a 
validated geopotential error model. To assess the error of the field by 
wavenumber, over different geographic regions, or in terms of orbit 
error, a detailed knowledge of model errors is required. Independent 
data needs to be obtained and calibration techniques need to be 
improved. This is especially important when assessing the separation of 
oceanographic and gravitational signals simultaneously contained within 
satellite altimeter observations. It is of critical importance when 
evaluating the accuracy of models for reaching the precision orbit 
determination (POD) performance objectives of future missions like 
TOPEX. In this latter case, actual observations are absent, and 

assessments rely solely on statistical considerations. Those 
statistics, which is found in the calibrated full solution covariances, 
need to be reliable. The uncertainty reported for the GEM-T1 

geopotential solution which are to be verified herein have been given as 
an RMS coefficient error by degree in Figure la and for each coefficient 
individually in Figure 1b. 

Obtaining detailed knowledge of the geopotential errors is made 
difficult due to the lack of a perfect force model, measurement model, 
and correction algorithms which enter into the solution. Although it 
may appear that we are limited to providing a single overall scaling 
factor (f in eq. 3-1) to account for systematic errors (a ) per data 
type (t) in order to obtain realistic error statistics, the factor 
o t may be varied somewhat to account for variations in f for different 
data types (see Appendix A2.6). We know that each error is unlikely to 
contaminate the model uniformly in a way representable by a single scale 
factor. However, only through a detailed and rather complete analysis 
directed at a total field calibration can the effects of each error 
source be isolated through the resulting behavior of subset 
calibrations. Appendix A3 derives an optimal weighting process where 
each data type is adjusted for a weighting factor f t . 
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In earlier gravitational modeling efforts, the inversion of the 
2500x2500 parameter matrix for GEM-T1 would have represented a 
considerable computational burden. On the IBM 360/95, which was the 
fastest mainframe available at GSFC until the late 1970' s, such a 
solution would have taken in excess of 100 CPU minutes. In 1982 the 
Cyber 205 vector computer was installed at GSFC. Our geopotential 
modeling efforts have benefitted substantially from the advent of this 
new "supercomputing" environment. On the Cyber, the inversion of the 
GEM-T1 normal matrix requires less than 2.5 minutes of CPU time. Given 
this enhanced capability, it is now possible to compute a variety of 
test solutions complete to degree and order 36 which will calibrate a 
broad spectrum of the errors in GEM-T1 . If this extensive testing 
yields comparable results to those which we found earlier (Marsh et al., 
1988), then we will have a more reliable verification of our estimates 
of the accuracy of the GEM-T1 model. 

In Table 3 the major data subsets of the GEM-T1 solution, as well 
as independent data used for calibration, are reviewed. The GEM-T1 
normal matrix is formed by summing the modularized contribution of each 
of these observation groupings. It is therefore a relatively simple 
task to manipulate these contributions and solve for models which either 
lack specific subset contributions or combine GEM-T1 with new 
independent observations such as altimetry and surface gravimetry. 
GEM-T2 and GEM-T3 are under development and will be models which 
optimally combine GEM-T1 with new data in this way. 

From equations (4.14) through (4.18) we have defined two methods 
for calibrating model errors from potential coefficient solutions; the 
first when the fields F and F are independent and share no common data; 
the second when the data within F are wholly contained within F. The 
test solutions made in conformance with the requirements of these 
calibration methods are given in Table 4. 

For each individual spherical harmonic, the calibration factors 
(eq. 4.18) tend to have random variations. To better assess the overall 
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TABLE 3 


MAJOR DATA SETS USED 
FOR THE CALIBRATION OF GEM-T1 
COEFFICIENT ERRORS 


MAJOR DATA SUBSETS WITHIN THE GEH-T 1 SOLUTION 
• SATELLITE TRACKING DATA 
LASER - 10 SATELLITES 
DOPPLER - 2 SATELLITES 


OPTICAL - II SATELLITES 


• COLLOCATION CONSTRAINTS - 
MODIFIED KAULA'S RULE 


r 10-5 

L /2 X 2 


] 


INDEPENDENT DATA COMBINED WITH GEM-T I FOR 
CALIBRATION PURPOSES 


• SEASAT ALTIMETER 

GLOBAL COVERAGE (OCEANS) 

• SURFACE GRAVITY 

OSU '86 NORMALS 
GLOBAL COVERAGE 

a 60 DEEP RESONANT SATELLITE 
LUMPED HARMONICS 
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Table 4 


Independent and Dependent 
Subset Solutions* 
for GEM-T1 Error Calibration 


Independent Data Solutions 

• GEM-T1 vs. Surface Gravimetry + SEASAT Altimetry 

• GEM-T1 minus LAGEOS vs. Surface Gravimetry + 
SEASAT Altimetry + LAGEOS 

Subset Data Solutions 


GEM-T1 vs. GEM-T1 minus 

• The laser data from GEOS-1, GEOS-2, GEOS-3 andBE-C; 
this data set is referred to elsewhere as "4-LASERS". 

• The laser data from STARLETTE. 

• The SEASAT and OSCAR Doppler data. 

• The optical observations acquired on 1 1 different satellites. 

• The laser data from LAGEOS 

GEM-T1 as a Subset Solution 

• Combining GEM-TI's data with surface gravimetry. 

• Combining GEM-TI's data with surface gravimetry and SEASAT 
alitmetry. 

• Combining GEM-TTs data with satellite deep resonant spherical 
harmonic constraints. 


* All solutions are complete to degree and order 36x36 and solve for 
ocean tidal and Earth orientation parameters. 
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quality of the calibration, we present the average calibration factors 
for terms of the same degree and terms of the same order separately. 
These are the k & and values given in eqs. (4.16) and (4.17) 
respectively. The size of the sample used to calculate these average 
calibration factors contributes to their overall stability whether they 
are collected by degree or order. Those sampling a larger set of terms 
tend to be better behaved across the spectrum within the model. This is 
probably why more stable results are generally found in the assessment 
by order for the low order terms and by degree for the high degree 
terms. Anomalous behavior of the calibration factors for averages taken 
by order are particularly useful for locating aliasing problems within 
the gravitational model arising from satellite dynamic sensitivities to 
terms above the 36x36 limits of the field. Performance of data sets 
from satellites experiencing larger atmospheric drag and other non- 
conservative forcing effects are also scrutinized to assess possible 
field contamination arising from these error sources. 

Briefly, Subsection 5.1 details the calibration of the GEM-T1 
solution, specifically the optimization of data weights, through a 
calibration of the solution variances. Gravitational models formed from 
completely independent data are evaluated by comparing their respective 
uncertainties with their coefficient differences. Subsection 5.2 shows 
a more direct calibration of GEM-T1 through a mutual calibration 
analysis where GEM-T1 is calibrated against other models which either 
lacked specific data subsets or added additional data to the GEM-T1 
field. Subsection 5.3 tests the sensitivity of the calibration 
methodology to the data weighting factor, f given in (eq. 3.1). A field 
which falsely reports superior results by having f increased by a factor 
of ten is evaluated in this subsection. Section 5.4 provides an overall 
summary of this calibration technique and the results obtained when 
evaluating GEM-T1. 
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5.1 Gravitational Model Calibrations Using Independent Data 

Turning first to calibrations using independent data, Figure 3 
shows the calibration factors, k & and k m , obtained when comparing GEM-T1 
with a field computed solely from a combination of SEASAT altimetry 
(where the altimeter ranges are treated as tracking observations) and 
block areal mean values of surface gravimetry. The gravimetry were 
provided to us by Pavlis (1988) in the form of compatible normal 
equations. These normal equations from gravimetry were developed under 
the direction of R. Rapp of The Ohio State University and have been 
scaled for consistency with the GEM-L2 geopotential. Here they will be 
used to calibrate GEM-T 1 . As noted in Marsh et al. (1988) the data in 
GEM-T1 is largely independent of GEM-L2 . The calibration factors, as 
shown in Figure 3, cluster around the ideal where k^=1 and 1^=1. This 
indicates that the overall uncertainties estimated within these 
independent models agree well with their coefficient variations. 
However, one of the limitations of this particular calibration test is 
revealed in Figure 4, which shows a large disparity in low degree 
coefficients between GEM-T 1 and this model computed solely from 
altimetry and surface gravimetry. A better test of GEM-TI's long wave- 
length harmonics was made by removing the Lageos data from GEM-T 1 and 
combining them with the normals of the altimetry /gravimetry model. The 
resulting two models were again composed of independent data, but both 
were now more equal in accuracy at low degree and closer to the uncer- 
tainty level of GEM-T1 itself. In Figure 5 the coefficient differences 
between GEM-T 1 lacking Lageos and this altimetry/gravimetry/Lageos 
solution are compared. The two solutions resulting from this shift of 
the LAGEOS information are compared in terms of their estimated 
uncertainty and rms coefficient differences in Figure 6. It is 
interesting to note that these uncertainties reverse if the optical data 
were also removed from GEM-T 1 and included in the independent test 
field. 


These "Lageos-shifted" models were then calibrated with one 
another. Figure 7 gives the average calibration factors (k^ and 1^) 
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which were obtained. Again, the results are well centered about 1, 
indicating that the model uncertainties are reasonably well 
calibrated. Furthermore, the calibration factors are more stable than 
those shown in Figure 3 although only marginally so. In summary, for 
these calibration procedures where geopotential models were composed of 
independent data sources, excellent overall calibration results have 
been obtained. The variances and model differences are consistent. These 
results indicate that the error estimates shown for GEM-T1 in Figures 
1.1 and 1.2 are supported by independent global altimetry and surface 
gravimetry . 

The large variation of from unity for high orders around m=28 
is due to the truncation at degree 36 for satellite resonant 
perturbations. Hence, this variation will consistently show up in the 
calibration plots for and will be discussed later in more detail. 


5.2 Gravitational Model Calibrations Using Dependent Solution Subsets 

Returning to Table 3, tracking data from 17 different satellites 
contribute data to GEM-T1 . These observations have been divided into 
several independent subsets. For each experiment which follows, the 
data was divided in such a way that the subset contributed significantly 
to the GEM-T1 solution. In some cases this required combining observa- 
tions from several satellites into a single subset. A solution was made 
after removing each of the data subsets individually, and all of these 
subset fields were used in an extensive evaluation of the uncertainties 
published for the GEM-T1 coefficients. 

A typical example of the variations found among the individual 
harmonic calibration factors, k^, is shown in Figure 8 for the case 
where GEM-T1 is calibrated against the model which removed the laser 
data contributions from GEOS-1,-2,-3 and BE-C. As expected, there is a 
good deal of variation in these calibration factors when going from 
term to term. However, the values overall fall into an approximate 
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FIGURE 6 



range of 0.2 to 2.0 and seldom depart by more than a factor of two from 
the degree and order average values which are shown in Figure 9. With 
the exception of the high order resonances of order 26 through 29, these 
average calibration factors are well behaved and centered about 1. As 
mentioned previously, the misbehavior at the high orders indicate that 
these satellites are sensitive to resonance terms beyond the 36th degree 
at these specific orders, and that there is some aliasing from these 
unmodeled higher degree contributions. When altimetry and surface 
gravity are introduced into the solution, these satellite data sets will 
contribute to the determination of higher degree resonance terms out to 
degree 50 and this source of aliasing will be eliminated. The normal 
equations from these and the other laser satellites are complete to 
degree 50 for all zonal, 1st and 2nd order, and resonance terms although 
they have not been used in GEM-T1 beyond degree 36. 

These calibration tests have been repeated for each of the subset 
solutions. Table 5 summarizes the overall calibration factors (defined 
as the mean of the k^ and k^ terms) from these assessments. With the 
exception of the model leaving out Lageos contributions, the overall 
calibration factors were all found to be highly consistent and close to 
the ideal value of 1. The less consistent Lageos results are somewhat 
understandable. The Lageos satellite orbits at nearly 6000 km altitude 
and senses little of the gravitational field beyond degree 8. It is 
easy to understand how this subset solution may be ineffective for 
calibrating a full 36x36 model. The low degree portion of the field (to 
6x6) however, gives essentially the same Lageos calibration factor of 
1.2 which is satisfactory. Also of note, the Lageos calibration is 
sensitive to the polar motion adjustment, and the calibration improves 
if the pole position is fixed. Also in Figure 7 the solution for Lageos 
with altimetry and surface gravity data result in a consistent 
calibration with the solution for GEM-T1 minus Lageos data. 

Additional tests were performed where GEM-TI's data became a 
subset of larger solutions— combination gravitational models which are 
being investigated as preliminary versions of GEM-T2 and GEM-T3. At 
this time, these more comprehensive models were solved only to degree 
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TABLE 5 


SUMMARY OF RESULTS 
1 FOR ERROR CALIBRATION 

CALIBRATION 

• GEM-TI vs GEM-TI minus DATA S1IRSFT FACTOR 


4 LASER SATELLITES (GEOS 1,2,3, BE-C) I 06 

STARLETTE LASER , 10 

OSCAR ♦ SEASAT DOPPLER , og 

OPTICAL (11 SATELLITES) 0 84 

LAGEOS LASER , 45 

• GEM-TI vs GEM-TI ♦ SURFACE GRAVITY 0 95 

• GEM-TI vs GEM-TI ♦ SURFACE GRAVITY + 

SEASAT ALTIMETRY 0 94 

• GEM-TI vs GEM-TI ♦ 

SATELLITE RESONANCE DATA 105 

• GEM-TI vs. SURFACE GRAVITY ♦ SEASAT 

ALTIMETRY 0 99 

• GEM-TI minus LAGEOS vs. LAGEOS + 

SURFACE GRAVITY + SEASAT ALTIMETRY 0 95 
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and order 36. Before the TOPEX work is completed, the limits of these 
fields will be extended to degree and order 50. These combination 
models were developed by adding new normal equations which were 
produced from SEASAT satellite altimetry and surface gravimetry to the 
GEM-T1 solution and truncating the final solution through parameter 
suppression at degree 36. Figure 10 presents one such calibration 
comparing GEM-T1 with the GEM-T1 /surface gravimetry combination 
solution. Again, a good calibration centering about 1 was obtained. 

Another source of independent satellite information was also used 
to calibrate GEM-T1 . Lumped harmonics representing the long term orbit 
changes of satellites undergoing deep resonant passages were utilized. 
In the first calibration method we compared the observed lumped 
harmonics directly with computed values of them using GEM-T1 . All of 
these observations were from unique orbits not found in GEM-T1. The 
observations used, their estimated accuracies and their projected GEM-T1 
errors are listed in Table 6. These observations provide a reliable 
error estimate for GEM-T1 since only lumped harmonics were used with 
errors considerably less than projected for GEM-T1 . We have computed an 
average calibration factor for GEM-T1 errors necessary to resolve the 
residuals of these observations with GEM-T1 predictions - namely: 


K 


60 


l 

i = 1 

(0.-C.) 

i l 


60 

v 


L 

i = 1 




1/2 


where Ch, aCL, (\ and oC. are respectively the observations, their 

estimated errors, GEM-T1 computed quantities and GEM-T1 predicted errors 

based on its variance-covariance matrix. AT. is the estimated 

truncation error for the ith observation for a model truncated at degree 

1 

36. Over this set of observations (note that aO.^ _< ^ aC^* for all i) we 
find that 
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TABLE 6. COMPARISON OF GEM-T1 WITH THE BEST INDEPENDENT SATELLITE 
RESONANT HARMONICS FOR CALIBRATION PURPOSES 






(1) 

(2) (3) 

X 

10 9 
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(6) 

(7) (8) 




RESONANCE 
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GEM-T1 







REF 










Satellite 



m 

q 

k 

Obs 

o Comp. 


AT 

RES 

RES/ RES/ 

Name 

Incl . 

Ecc. 







c 


(O-C) 

°c e 




(1) 

11 

-2 

1 

C -38.1 

1.0 -45.7 

5.6 

1.0 

7.6 

7.6 1.28 

VANGUARD 3 

33.35 

.188 





S -33.4 

1.2 -17.6 

5.6 

1.0 

-15.8 

-13.2 -2.64 

VANGUARD 3 


(1) 

11 

-1 

1 

C -30.6 

1.0 -35.0 

5.2 

0.6 

4.4 

4.4 0.80 

VANGUARD 3/2R 

33.13 

.184 





S 0.3 

2.1 1.9 

5.3 

0.6 

-1.6 

-0.8 -0.27 

VANGUARD 3/2R 




11 

0 

1 

C 26.1 

2.2 25.9 

6.6 

0.6 

0.2 

0.1 0.03 

VANGUARD 3/2R 







S 52.6 

1.1 53.7 

6.7 

0.6 

-1.1 

-1.0 -0.15 

VANGUARD 3/2R 




11 

1 

1 

C -42.9 

1.5 -40.7 

8.6 

0.8 

-2.2 

-1.4 -0.23 

VANGUARD 3/2R 







S -13.9 

2.9 -8.2 

8.0 

0.8 

-5.7 

-2.0 -0.64 

VANGUARD 3/2R 




11 

2 

1 

C 13.6 

1.5 6.5 

9.9 

1.2 

7.1 

4.8 0.68 

VANGUARD 3/2R 







S 24.5 

1.4 47.0 

10.8 

1.2 

-22.5 

-16.1 -1.96 

VANGUARD 3/2R 



(1) 

12 

0 

1 

C -4.2 

2.0 -0.5 

7.0 

0.9 

-3.7 

-1.7 -0.48 

D- 1 A & D-1 A ROC 

34.10 

.139 





S 38.4 

1.4 40.4 

7.0 

0.9 

-2.0 

-1.2 -0.27 

D- 1 A & D-1 A ROC 




12 

1 

1 

C -42.4 

1.4 -29.2 

10.3 

1.2 

-13.2 

-7.2 -1.20 

D-1A & D-1 A ROC 







S -23.4 

2.4 -15.3 

9.9 

1.2 

-8.1 

-3.0 -0.76 

D-1A & D-1 A ROC 



(1) 

12 

1 

1 

C 19.2 

0.5 16.9 

4.1 

0.5 

2.3 

3.1 0.52 

EX. 9 ROC. 

38.90 

. 121 





S 42.2 

1.5 39.1 

4.2 

0.5 

3.1 

2.0 0.67 

EX. 9 ROC. 



(1) 

12 

0 

1 

C -44.3 

2.1 -34.1 

6.2 

0. 1 

-10.2 

-4.8 -1.48 

TIROS 9 

96.40 

.117 





S 101.9 

1.3 100.8 

6.3 

0.1 

1.1 

0.8 0.16 

TIROS 9 



12 

1 

1 

C -1.1 

1.2 2.3 

3.9 

0.1 

-3.4 

-2.8 -0.80 

TIROS 9 







S -14.0 

1.1 -23.1 

3.8 

0. 1 

9.1 

8.2 2.20 

TIROS 9 



(1) 

14 

0 

1 

C -31.6 

1.1 -31.0 

9.5 

1.3 

-0.6 

-0.4 -0.06 

0V2-1 ROC. 

32.10 

.013 





S -4.8 

1.1 17.8 

9.8 

1.3 

-22.6 

-13.2 -2.17 

FRAGMENTS 



(2) 

14 

0 

1 

C -1.8 

0.6 -2.0 

2.2 

0.3 

0.2 

0.2 0.07 

METEOR 10 

81 .23 

.007 





S -19.5 

0.5 -17.0 

2.2 

0.3 

-2.5 

-4.3 -1.06 

METEOR 10 



(1) 

14 

0 

1 

C -2.3 

0.3 -0.3 

2.3 

0. 1 

-2.0 

-6.2 -0.83 

LANDSAT 1 

99.10 

.001 





S -22.5 

0.5 -24.4 

2.3 

0.1 

1.9 

3.8 0.80 

LANDSAT 1 



(2) 

14 

0 

1 

C -2.8 

0.3 -0.3 

2.2 

0.1 

-2.5 

-7.8 -1.05 

LANDSAT 1 

99.10 

.002 





S -22.6 

0.1 -24.4 

2.2 

0.1 

1.8 

12.1 0.77 

LANDSAT 1 



(3) 

15 

0 

1 

C -26.3 

1.7 -3.1 

13.0 

3.0 

-23.2 

-6.8 -1.65 

1 965-9 A 

31.76 

.007 





S -11.5 

0.8 4.5 

13.0 

3.0 

-16.0 

-5.1 -1.14 

1 965-9 A 



(3) 

15 

0 

1 

C 26.3 

1.0 2.4 

12.4 

2.9 

23.8 

7.9 1.78 

1 969— 68b 

32.97 

.004 





S 8.1 

1.2 -4.9 

12.4 

2.9 

13.0 

4.2 0.97 

1969-68B 



(1) 

15 

0 

1 

C 28.3 

1.5 2.4 

12.2 

2.9 

25.9 

8.0 1.96 

TETR-3 

33.10 

.012 





S 7.4 

1.5 -4.9 

12.3 

2.9 

12.3 

3.8 0.93 

TETR-3 



(3) 

15 

-1 

1 

C -46.5 

2.7 -50.1 

11.2 

1.8 

3.6 

1.1 0.30 

1970-1 1 1 A 

74.00 

.001 





S -40.5 

4.0 -46.5 

11.2 

1.8 

6.0 

1.4 0.48 

1970-1 1 1 A 




15 

0 

1 

C -26.0 

1.0 -21.2 

6.0 

0.6 

-4.8 

-4.2 -0.75 

1970-1 1 1 A 







S -5.2 

1.3 -9.8 

6.0 

0.6 

4.6 

3.3 0.72 

1970-1 11 A 



(3) 

15 

-1 

1 

C -45.5 

2.0 -50.2 

11.3 

1.8 

4.7 

1.8 0.39 

197 1 — 1 3 1 B 

74.05 

.002 
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S -35.2 

1.0 -46.6 

11.3 

1.8 

11.4 

5.6 0.95 

1971-131B 



15 

0 1 C -24.6 

1.3 -21.3 

6.0 

0.6 

-3.3 

-2.3 -0.51 

1971 - 1 3 1 B 




S -6. 1 

1.0 -9.8 

6.0 

0.6 

3.7 

3.2 0.58 

1971-131B 



(3) 15 

0 1 C -23.0 

1.6 -20.7 

5.5 

0.3 

-2.4 

-1.5 -0.39 

1967-42A 

80.17 

.007 


S -8. 6 

1.3 -6.9 

5.5 

0.3 

-1.7 

-1.3 -0.29 

1967-42A 



15 

1 1 C -54.6 

3.2 -45.6 

10.8 

0.9 

-9.0 

-2.7 -0.76 

1967-42A 




S -37.1 

2.6 -45.7 

10.7 

0.9 

8.6 

3.1 0.74 

1967-42A 



(3) 15 

-1 1 C -62.9 

2.6 -49.0 

13.7 

1.2 

-13.8 

-4.8 -0.95 

1971 —54a 

90.23 

.002 


S -53.4 

1.6 -56.3 

13.6 

1.2 

2.9 

1.5 0.20 

1971 -54a 



15 

0 1 C -16.0 

0.2 -15.8 

2.7 

0. 1 

-0.2 

-0.8 -0.07 

1971 —54A 

90.21 

.002 


S -6.9 

0.2 -6.8 

2.6 

0.1 

-0.1 

-0.5 -0.04 

1971 -54a 



(4) 27 

0 1 C -14.5 

2.0 -13.4 

6.2 

1.4 

-1.1 

-0.5 -0.16 

NAVSATS 

89.80 

.010 


S +10.0 

2.0 6.8 

6.1 

1.4 

3.2 

1.3 0.47 

(1967-48A, 










70-67A, 73- 

-81 A) 


(5) 28 

0 2 C 7.0 

1.1 2.9 

4.7 

1.0 

4.1 

2.8 0.80 

METEOR 10 

81 .23 

.007 


S 1.5 

1.1 0.4 

4.7 

1.0 

1.1 

0.7 0.21 

METEOR 10 



(6) 28 

0 2 C -9.3 

0.9 -14.3 

9.2 

1.6 

5.0 

2.7 0.51 

LANDSAT 1 

99.10 

.002 


S -21.7 

0.9 -7.4 

9.2 

1.6 

-14.3 

-7.8 -1.45 

LANDSAT 1 



(3) 30 

0 2 C -9.8 

0.6 0.6 

8.2 

0.4 

-10.4 - 

14.8 -1.21 

1971 -54A 

90.23 

.002 


S 9.0 

0.7 0.1 

8.2 

0.4 

8.9 

10.5 1.03 

1971 —54A 



(7) 31 

0 2 C -2.9 

1.2 0.4 

6.9 

2.2 

-3.3 

-1.3 -0.43 

1961-o1 

97.20 

.005 


S 9.0 

2.2 -2.4 

6.9 

2.2 

11.4 

3.6 1.44 

1961-qI 



Column 

Notes : 









( 1 ) i (2) 

For all of 

these observation 

pairs 

CVS: o 

0 (C) . o 0 (S> 

< (o c (C) + 

o s (C)]/3. 


Thus the comparisons are 

all 

dominated by 

the estimated model 



hence they serve principally to 

calibrate the model. 


error o^, and 



fThe resonances are all of periods > 14 days with most passing through perfect 
commensurability with the earth’s rotation. Their frequencies are: 

^m,q,k = wq ♦ (oh-M) k + m (fl-e e ) 

where msgeopotential order and M,n,w are the satellite’s mean motions of mean 
anomaly, argument of perigee and node. e e is the earth's rotation rate. C and S are 
the cosine and sine terms of the lumped harmonic observed for these resonances. The 
definition of these harmonics is from Ref. (1). Harmonics from all other references 
have been converted to this definition which normali 2 es the observation with respect 
to the geopotential harmonic (C ) of greatest sensitivity to the orbit element whose 
resonance is observed. Here the elements are mainly inclination but equivalently 
semi-major axis eccentricity and along track perturbation, covering all cases. 

(3), (4) Computed observation and estimated commission error from GEM-T1 and its 
covariance matrix as presented in this report (precalibrated or formal errors) 

(5) Estimate of truncation error for resonant harmonic of GEM-T1 for 
terms fc>36 using a power law for the geopotential signal . 


D - (3) _ |1/2 
2 


/o 


tr 


(6) Residual = 

(7) Residual * 
for the terms in _its solution 
(6) Residual * 


* AT 2 ♦ {K o c ) 2 


: weighted residual if model were perfect 
1/2 


: weighted residual accounting 


for calibrated commission error. 

Here the calibration factor was found by solving for K from 


60 


7 REST * (a . + Lt. 
i=1 1 °' i 1 


K 2 o l^) = 60 


The result was : K=1.05. The resulting statistics in (8) 
yield : #>2o : 3 (expected from unit normal : 3) 

#>1o : 17 (expected from unit normal : 19) 
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K = 1.14. 


To test this calibration of the lumped harmonics of GEM-T1 with 
the dependent data set method in terms of the specific geopotential 
harmonics in GEM-T1 (instead of the lumped harmonic observables) we have 
combined the lumped harmonics data with GEM-T1 . This new solution was 
then compared to GEM-T1. As expected from the small set of 60 lumped 
coefficients, only a limited number (about 280) of resonant harmonics 
(fi,,m) in GEM-T1 were changed appreciably by these new data. The rms 
average calibration factor K (Eq. 4.18) for these 60 observations was 
found to be 1.05 and agrees well with the 1.14 calibration factor 
obtained above for the 60 lumped harmonics. The importance of these 
resonance tests is that they show the capability of GEM-T1 over a wide 
range of "new orbits" with respect to the dominant perturbations of 
those orbits. Furthermore, since the computed errors for these lumped 
harmonic "observations" used the full variance-covariance matrix of 
GEM-T1 , these are direct tests of the validity of that matrix, a task 
essential to reliable radial error prediction for TOPEX. Many more such 
full matrix calibrations will be presented in Section VI. 

Table 5 (last five entries) also summarizes the overall 
calibration factors which were determined for GEM-T1 from these tests 
using more advanced gravitational models. In all cases tested, the 
overall calibration factors indicated satisfactory error estimates for 
both GEM-T1 and these more advanced gravitational field coefficients. 

5.3 Tests of Overall Observation Scaling Factor 

Tests were also made to determine the sensitivity of our 
calibrations to the overall scaling factor f of the solution (eq. 
3.1). Clearly, the overall formal uncertainties of the model are 
strongly affected by the adopted value of f used in combining the signal 
and observation matrices. By design in GEM-T1 and previously in GEM-9 
(Lerch et al., 1979), the scaling of f was found to produce solution 
uncertainties which reflected the true accuracy of the field. For 
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GEM-T1 we found f=.02. To test this assignment using our more 
comprehensive calibration methods both GEM-T1 and the GEM-T1 subset 
model which lacked data from the 4-Laser satellite group were recomputed 
with f increased by a factor of 10 to f=0.2 while leaving the signal 
matrix weight for the Kaula constraints (o ) unchanged. This increased 

A / 

the satellite data weight in the solution by a factor of 10, and as 
expected, the overall parameter standard deviations improved by 
approximately /10. This was expected since, as we previously mentioned, 
(Section III) tests with global sets of gravimetry and other independent 
data gave f=0.02 as the optimum scale for the satellite data relative to 
the a priori constraint information. Here, we deliberately corrupted 
this weighting and attempted to confirm this estimate with reference 
only to the solution's own satellite contribution. 

As shown in Figure 11, when calibrated using the subset solution 
technique, the calibration factors k & and ^ were both found to yield 
values approximately equal to /10 for these test models indicating that 
*02 is indeed close to optimum. The subset solution method is thereby 
shown to be capable of revealing both a poorly calibrated model and also 
showing the correct approximate level at which these errors were 
improperly reported. 

5.4 Summary of the Calibration Based Upon GEM-T1 Variances 

Clearly from the results shown within this section, the calibra- 
tion method we have developed is sensitive to overly optimistic solution 
statistics and is consistent across many data subsets as they are 
individually tested against the complete GEM-T1 field. These results 
therefore confirm the coefficient uncertainties of GEM-T1 for the 
diagonal terms of the covariance matrix. We will now proceed to test 
the full covariance matrix for this solution. 
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GEM-T1 with 10 times the DATA WEIGHT vs. GEM-T1 minus 4-LASER with 10 times the DATA WEIGHT 
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DEGREE / ORDER 

CALIBRATION OF ERROR ESTIMATES BASED UPON SUBSET SOLUTIONS 
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VI. GRAVITY MODEL CALIBRATION BY EIGENVALUE-EIGENVECTOR ANALYSIS 
6.0 Introduction 

One of the major objectives of this research effort is to produce 
gravitational models capable of supporting the stringent orbit 
determination requirements of the TOPEX Mission. However, to fully 
assess the performance of a gravitational model on a satellite orbit 
(prior to launch), knowledge of only the error variances of the field is 
insufficient to make a reliable projection. For example, consecutive 
even and odd zonal harmonics are correlated by .95 and significant 
correlation exists for coefficients in all orders (particularly 
resonance orders). For TOPEX, the predicted radial error based upon the 
GEM-T1 covariance matrix (as more fully described later) gives about 
25 cm while using the variances alone predicts radial errors for TOPEX 
exceeding 80 cm RMS. 


Any given satellite orbit samples the earth's gravitational field 
in a way which causes it to sense certain perturbative frequencies. 
Each of these perturbations can be mathematically described as some 
linear combination of the spherical harmonics used to represent the 
gravitational field (see for example Wagner and Klosko, (1977) or 
Reigber, (1974) who base their analyses on the linear perturbation 
theory given in Kaula, (1966)). These sums (or "lumped-harmonics") can 
be very accurately determined, even more so than the individual harmonic 
coefficients, but they are satellite-specific and tend to reflect the 
distribution of orbital characteristics sampled when forming the 
model. Therefore, our pre-launch analysis was directed at estimating 
the errors in the "lumped-harmonics" which are TOPEX specific and 
TOPEX sensitive. In particular, we are most interested in those that 
represent significant radial perturbations. 

WbiiG it is correct to assume that highly correlated parts of our 
gravitational model have degraded variances due to significant off- 
diagonal terms in the normal matrix, only these cross-correlations allow 
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one to properly sum the effects within a given order to calculate errors 
in these "lumped-harmonics'*. Specific knowledge of the full covariance 
(not merely the standard errors) of the solution, especially the cross- 
correlation of the coefficients with odd and even parity within the same 
orders is therefore essential to predict orbit errors. Acknowledging 
this, we wondered, for example, whether the Kaula constraint matrix 
within our solution (which artificially reduces the cross-correlation 
between high degree terms) would yield a seriously distorted covariance 
matrix. In any case, it is necessary to provide a calibration of the 
complete GEM-T1 covariances before covariant orbital accuracy 
projections on TOPEX can be considered reliable. 

Our approach to the covariance calibration of GEM-T1 was through 
the eigenvectors of its error covariance matrix. Again, we were 
fortunate that we had a supercomputer available for calculating the 
eigenvectors within the GEM-T1 solution. The full covariance matrix for 
each of the solutions described in Table 4 were calculated and provided 
the means for calibrating the GEM-T1 covariance. First, in Section 6.1, 
we review the mathematical development of this calibration technique 
which strongly parallels that given in Section IV. We then discuss the 
projection of the eigenvectors from various dependent models onto those 
of GEM-T1 using the differences in their respective covariances (Section 
6.2). Section 6.3 shows the calibration of the model using all of the 
major data subsets within GEM-T1 and results obtained after introducing 
surface gravimetry and satellite altimetry into the model. Here an 
approach which projects the difference of the coefficients directly on 
the eigenvectors of the main or subset field is employed. 

6.1 Calibration by Eigenvector Analysis: Direct Method 

Given two fields F and F, let their corresponding coefficient 
solutions be C and C, and their corresponding error-covariance matrices 
be V and V, respectively. 
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Since we are interested here in assessing only the gravity model 
uncertainties and their effects on orbital trajectories, the full error- 
covariance matrices were truncated to include only the geopotential 
parameters. Thus the eigenvector technique of calibration (to be 
discussed herein) considers the previously ignored covariances. 

Let x and X be the diagonal eigenvalue matrices of the now 
truncated V and V matrices, respectively. By convention, x and \ are 
arranged in descending order. Thus, the matrices take on the form 


*1,1 0 o 

0 * 2,2 0 

0 0 *3,3 

0 0 0 

0 0 0 


0 

0 

0 

*4,4 

0 


0 

0 

0 

0 


••• *1363,1363 I 


( 6 . 1 ) 


Let B and B correspond to the ordered eigenvector matrices, where 
the individual normalized eigenvectors form the columns of B and B. 
Therefore, 


V = B X B T 


(6.1a) 


V : B X B T 


(6.1b) 


The eigenvectors of the error-covariance matrix represent the 
independent combined harmonics unique to the solution while the 
corresponding eigenvalues represent the variance of the error in 
estimating each of these combinations. These combined harmonics C' are 
defined from C as: C' = B^C. Thus, a small eigenvalue means that a 
particular combination (corresponding approximately to a "lumped" 
harmonic of all terms in the field) is well-determined by the solution 
while a large eigenvalue means that a particular combination is poorly- 
determined . 
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Now we can derive eigenvector calibration factors in analogy to 
equation (4.1) 


fi,,m 


= AC / 

fl,,m 


'% , m 


Let 


AC = C 


C = ( 5 - ' f C - C true) 


( 6 . 2 ) 


A V = V - V 


(6.3a) 


when the data included in F is a subset of the data included in F. 
Also, again by analogy, when models are compared based exclusively on 
independent data, 

AV = V + V (6.3b) 


By equation (6.1), where the D subscript represents the difference in 
the covariances as in (6.3a), we get 


AV 


b d x d b d 


( 6 . 1 )) 


While equation (6.4) is similar to equation (6.1), the inter- 
pretations of the eigenvalue sizes are quite different. Small eigen- 
values derived from V or V indicate that the associated eigenvectors 
are well-determined by the data used to form the respective F or F 
gravity fields. In contrast, small eigenvalues derived from AV imply 
that the parameters associated with the eigenvectors are not signifi- 
cantly altered due to the presence or absence of a specific data subset 
which delineates the difference between the main (F) and subset field 
solutions (F). Hence, significant effects of the data subsets are found 
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with the largest eigenvalues of X D . In this case, the magnitude of the 
eigenvalues is a measure of the improvement in the variances achieved 
through the inclusion of these data. 

The direct eigenvector projection Ay of the difference in the 
actual harmonics AC of the two solutions can be obtained through the 
eigenvector matrix by the relation 

Ay = AC (6.5) 

We shall refer to Ay as the eigenvector parameters whose variance is 
T 

given by E(Ay Ay ) : 


EUy iy T ) = Bp E( AC AC T ) B Q (6.6) 

We are now in a position, by analogy to (4.1), to define the 

eigenvector calibration factors where AC becomes I Ay I while e 

ip ^ /p x , m am 

becomes (E(Ay Ay )) as follows: 

CF. = | Ay | / [ E ( Ay Ay ) T ] 1 72 


From (4.3) and (4.4) we have shown that 
E( AC AC T ) = AV 


Thus , 


E(Ay Ay T ) = AV B D 


Now 


b d b d = b d B l = 1 


(6.7) 


( 6 . 8 ) 


(6.9) 
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from the definition of the eigenvector matrix. Using these identities, 
equation (6.4) can be transformed to 


aV B D = *D 


( 6 . 10 ) 


so that 


E(Ay Ay T ) = X^ 


( 6 . 11 ) 


The direct eigenvector calibration factor then becomes 

CF. = [(Ay Ay T ).. / (X D )..] 1/2 (6.12) 

where i = 1 to 1363, the number of geopotential parameters in GEM-T1 
type solutions. 

Table 7A shows the results of three calibration tests. For 
simplification purposes and to average over individual samples, the 
individual calibration factors are divided into groups of 50 beginning 
with the CF^ associated with the largest eigenvalue. An RMS of the 
CF.'s 


<CF. > 
k 


n+49 

I 


i= n 


(CF.) 2 / 50 


1/2 


(6.13) 


and an average of the <X n > 

D k 


<V. 


n+49 

I 

i=n 


< X D ) ii 7 50 


(6.14) 


where n = 50(k-1) + 1, are obtained for each group. Finally, an rms of 
the CF k ’s (k = 1 to 27) for the complete sets of CF i 's, 
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rms 


1/2 


(6.15) 


- 27 


- k=1 


<CF. > £ 
k 


/ 27 


and a weighted rms 


Wrms = 



<CF, > 2 
k 



( 6 . 16 ) 


are computed for each of the three tests. Note, the weighted rms gives 
added weight to the eigenvector parameters which are more significantly 
affected by the data. 


The exponential fall off in average eigenvalue size shown in 
Table 7A led directly to the computation of the weighted rms calibration 
factor in equation (6.16). Using this formula, the overall weighted 
calibration factor now more closely reflects the significant differences 
in the errors of the two solutions. 


The size of the average eigenvalues in Table 7A is a measure of 
how much a tested data set improves the main field over the subset model 
lacking these data. In general, the average eigenvalues for the surface 
gravity and altimeter data set are larger than the corresponding average 
eigenvalues for the other two tested data sets. Furthermore, the data 
subset from the 4 laser satellites ('’4-LASERS") improves the solution 
more than the solitary Starlette data set. Since four unique orbital 
geometries are sampled in the 4-LASER case, this result is expected. 

Note that there are only 25 (instead of 27) eigenvalue groups of 
50 samples for the case of [GEM-T1] vs. [GEM-T1 minus STARLETTE data]. 
The algorithm which computed the eigenvalues produced about 100 negative 
but small eigenvalues. Negative eigenvalues have no physical signifi- 
cance but reflect numerical properties of the computation environment. 

1 P 

The positive eigenvalues have a range of about 10 . The computer 
arithmetic, however, was done using only 14 digits. Therefore these 
negative eigenvalues undoubtedly result from machine round-off errors. 
Thus, for all practical purposes, they can be treated as zeros. The 
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1 ABLE 7. A 


DIRECT EIGENVECTOR CALIBRATION 
FACTORS <CF k > AND AVERAGE EIGENVALUES <X Dk > 


GEM-T1 vs. GEh-Tt Minus/Plus the following dote sets 


FIGENVALUE 

(-)STARLETTE 

GROUP NO(k> 

<C£ k > 

<Xj> k > 

1 

0.87 

.1 IE-15 

2 

1.15 

.60E-16 

3 

1.13 

.36E- 16 

4 

1.55 

.2 1 E- 1 6 

5 

1.24 

.95E- 17 

6 

1.00 

.33E- 17 

7 

0.72 

.10E-17 

8 

0.66 

.36E- 18 

9 

0.64 

.14E-18 

10 

0.52 

.53E-19 

1 1 

0.52 

.20E-19 

12 

0.44 

.B4E-20 

13 

0.45 

.36E-20 

14 

0.40 

.1 6E- 20 

15 

0.35 

.71E-21 

16 

0.38 

.31E-21 

17 

0.33 

.13E-21 

18 

0.39 

.54E-22 

19 

0.39 

.2 1 E- 22 

20 

0.27 

.74E-23 

21 

0.34 

.24E-23 

22 

0.27 

.76E-24 

23 

0.31 

.22E-24 

24 

0.28 

.49E-25 

25 

0.36 

.63E-26 


?6 

?7 


(-) 4 LASERS 


<CF k > 

<Aj>|c> 

0.92 

.1 1 E- 1 5 

0.66 

.54E-16 

1.10 

.36E-16 

1.26 

.24E-16 

1.03 

.1 5E- 1 6 

1.19 

.94E-17 

0.98 

.55E-17 

1.26 

.3QE-17 

1.08 

.1 5E- 1 7 

0.96 

.8 1 E- 1 8 

0.88 

.43E-18 

0.65 

.23E-18 

0.80 

.1 2E- 1 8 

0.60 

.66E-19 

0.82 

.35E-19 

0.60 

.19E-19 

0.59 

.10E-19 

0.73 

.52E-20 

0.58 

.27E-20 

0.64 

.1 4E-20 

0.59 

.67E-21 

0.58 

.31E-21 

0.63 

.1 3E-2 1 

0.56 

.51 E-22 

0.56 

.1 7E-22 

0.61 

.37E-23 

0.67 

.35E-24 


<♦> SURFACE GRAVITY 
♦ ALTIMETER 


<CF k > <Ax> k > 


0.73 

.28E- 15 

0.85 

. 1 4E- 1 5 

0.82 

.1 IE-15 

0.88 

.83E- 1 6 

0.95 

.68E- 16 

0.85 

.57E-16 

0.79 

.50E- 16 

0.83 

.45E-16 

1.16 

.40E- 16 

1.12 

.36E- 1 6 

1.06 

.32E- 16 

1.03 

.30E- 16 

0.84 

.27E-16 

0.93 

.24E-16 

1.04 

.20E- 16 

1.04 

. 1 6E- 1 6 

1.1 1 

.1 IE-16 

1.27 

.76E-17 

1.52 

.46E- 17 

1.69 

.26E-17 

1.74 

.14E- 17 

1.31 

.65E-18 

1.64 

.26E- 18 

2.03 

.8 1 E- 1 9 

1.85 

.1 9E- 1 9 

1.91 

.25E-20 

2.01 

.14E-21 


RMS 0.69 

WEIGHTED 1.08 

RMS 


0.89 

1.00 


1.30 

0.88 
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weighted rms calibration factors for the three test cases presented in 
Table 7 A show that the data in GEM-T1 is approximately weighted 
correctly. Calibration factors less than 1 mean that the errors have 
been over-estimated. Thus, for the cases of [GEM-T1 ] vs. [GEM-T1 minus 
STARLETTE data] and of [GEM-T1 ] vs. [GEM-T1 minus 4 LASER data], the 
attributed errors are actually slightly too large, but only if all the 
eigenvector parameters are weighted the same. 

For the case of [GEM-T1 plus surface gravity + altimeter data] vs. 
[GEM-T1], the rms calibration factor indicates that at certain 
frequencies the attributed errors have been under-estimated. However, a 
close look at Table 7A shows that all of the large calibration factors 
are associated with small average eigenvalues. In this case, the small 
eigenvalues generally correspond to eigenvectors of primarily long- 
wavelength combination terms or zonal and resonance effects. Surface 
gravity and altimeter data cannot resolve these dynamically resolved 
terms nearly as well as the direct satellite tracking (see Figure 4). 
Therefore, it is not surprising that these eigenvector calibration 
factors for the best determined eigenvectors are all significantly 
larger than 1 using surface gravimetry and the radially-sensitive 
satellite altimetry. Hence, the weighted calibration factor (0.88) 
should give the more reliable estimate. 

The average geoid height error Ah ± implied by each eigenvalue is 


A h i = /X ii R e (6.17) 
where R 0 is the mean radius of the earth. 

Thus, the contribution of these long wavelength eigenvalue errors is 
insignificant compared to the contribution of the large shorter 
wavelength eigenvalue errors. Therefore, if the attributed errors (for 
the surface gravity/altimeter case) for these longest wavelength terms 
were increased to reduce their calibration factors to approximately 1, 
there would be little change in the total predicted orbit determination 
error produced by GEM-T1 . 
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Of the calibration factors associated with large average eigen- 
values, only a few of them are exceptionally high, such as the 1.55 in 
the case of [GEM-T1 ] vs. [GEM-T1 minus STARLETTE data]. A close look at 
the components of the associated eigenvectors shows that the largest 
components are higher order resonant and near-resonant terms, such as 
C(33,27) and S(33,28). The primary resonance order of STARLETTE is 14 
with strong secondary resonances with 27th and 28th order terms. 
STARLETTE senses terms beyond the degree 36 cutoff of GEM-T1 for these 
resonance orders. Therefore it appears that these large calibration 
factors are due to aliasing arising from these neglected coeffic- 
ients. Thus, if the gravity field model is expanded to degree and order 
50 as is planned when altimetry and surface gravimetry are fully 
utilized, STARLETTE will contribute significantly to the recovery of 
these higher degree harmonics. Therefore, this aliasing problem and the 
large calibration factors for these vectors would be expected to 
disappear. 

When the weighed rms calibration factors are computed for each of 
the three cases, the results are all approximately 1. Thus, as far as 
the significant errors in the GEM-T1 gravity field model are concerned, 
especially for the more significant larger errors, a satisfactory 
calibration is indicated. 


6.2 Calibration by Eigenvector Analysis: Projection Method 

Calibration can also be done by an indirect projection method 
using eigenvector analysis. Instead of converting AC into combined- 
harmonics using the eigenvectors of the difference of the two error- 
covariance matrices by equation (6.5), AC can be projected onto the 
combined-harmonics of either the main gravity field solution or the 
subset gravity field solution using either B or B, respectively. 

Ay p = B T AC (6.18a) 
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(6.18b) 


Ay p = B T AC 
Thus , we have 

E (Ay p ^p) = B T AV B (6.19a) 

E (Ay p Ayp) = b t «v b ( 6 . 19 b) 

The calibration factors of the projection are then 

PCF. = [ ( Ay p iyp) / (B T AV B)..] 1/2 (6.20a) 

PCF = l ( Ay „ Ay ) / (B T AV B) ] (6.20b) 

ii 

The eigenvectors of any GEM-T1 -type field completely span the 1363 
geopotential parameter space. Thus, mathematically, Ay and its 
associated eigenvector calibration factors can be computed using either 
B or B. Projecting AC and AV into the combined-harmonics of either the 
main or subset fields (rather than onto the harmonics of the difference 
field) yields additional insight for assessing the errors in the gravity 
field solutions. 

Table 7B shows the results of the projection method when applied 
to the three test cases shown in Table 7A. Individual calibration 
factors were calculated for each eigenvector and averaged in groups of 
50 as before. However, for the sake of brevity, only the rms and 
weighted rms of the projected eigenvector calibration factors are 
presented in the table. 

The results show again that the data in GEM-T1 is basically well- 
calibrated. When projections are made onto the main field, both the rms 
and weighted rms projected calibration factors are found to be approxi- 
mately 1 . 
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TABLE 7.B 


PROJECTED EIGENVECTOR 
CALIBRATION FACTORS <PCF fc > 

GEM-T1 vs. GEM-T1 Minus/Plus the following data sets 


(-) STARLETTE 


PROJECTED ONTO 


(-) 4 LASERS 


PROJECTED ONTO 


(♦)SURFACE GRAVITV 
♦SEASAT ALT IMETRy 

PROJECTED ONTO 



SUBSET 

FIELD 

(GEM-T1) 

MAIN 

FIELD 

SUBSET 

FIELD 

(GEM-T1) 

MAIN 

FIELD 

(GEM-T1) 

SUBSET 

FIELD 

MAIN 

FIELD 

RMS 

1.30 

1.03 

1.15 

.99 

1.22 

.87 

WTED. 

RMS 

1.09 

.99 

1.00 

.94 

.89 

.97 
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However, when projections are made onto the subset field, the rms 
projected calibration factors depart markedly from 1. This result 
should not be surprising since the subset field does not contain any of 
the data whose errors are being assessed. Thus, the eigenvectors of the 
subset field are not likely to align with the errors of the data being 
calibrated, resulting in less reliable projections. On the contrary, 
projections onto the main field which contain the data being calibrated 
are likely to align with the data errors yielding more reliable 
calibration results. 

A closer look at the projected calibration factors when using the 
subset solution for each group of 50 also shows that all of the high 
calibration factors are associated with small average eigenvalues. The 
weighted rms projected calibration factors which favor more significant 
model differences are all reasonably close to 1. Therefore, as far as 
the significant errors in the data are concerned, the GEM-T1 gravity 
field is still calibrated correctly even when the errors are projected 
onto the subset field. We thereby have another useful method for 
further assessing the errors in GEM-T1 through this projection method 
utilizing subset solutions. 


6.3 Calibration by Means of Eigenvector Projections onto GEM-T1 

A total of 1363 projected eigenvector calibration factors were 
computed for each of the subset test solutions. To reiterate, the 
PCFji's were then ordered so that the PCF-j term came from the best 
determined eigenvector found within the GEM-T1 solution. Likewise, 
PCF i 363 was from the poorest determined eigenvector. To more easily 
evaluate the behavior of these calibration factors (as was done for the 
calibration factors k^ and k m ) , the CPF.^ factors have been averaged over 
groups of 50 taken from this ordering. The first group contained the 
best determined eigenvectors, while the 27th group had the worst. 
Table 8 presents the results of these average projected eigenvector 
calibration factors from tests using: 
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Case (A) a GEM-T1 subset solution that lacks the contribution of 

GEOS 1, 2, 3 and BE-C laser data: "4 LASERS" and also 

Case (B) GEM-T1 as the subset of a GEM-T1 plus altimeter/surface 

gravimetry model. 

Shown in Table 8 are the average size of the GEM-T1 eigenvalues 
and the scaled average error (found through a scaling with R e ) of the 
eigenvalues in cm within each of the groupings. The group calibration 
factors for each of these test solutions are also given in Table 8. In 
general, the projected eigenvector calibration factors have stable 
behavior over these groupings, and as seen earlier with the coefficient 
uncertainties, these eigenvector calibration factors also cluster around 
1 indicating GEM-TI's covariance is well calibrated. 

All of the subset solutions were evaluated using this projected 
eigenvector assessment including those which used subsets of GEM-T1 and 
those where new data was added to the GEM-T1 normals. Table 9 compares 
the weighted rms calibration factor obtained as in equation (6.16) from 
the eigenvector analysis with that obtained earlier using the coeffic- 
ient differences and their variances. 

The eigenvector analysis was also performed on the GEM-T1 plus 
deep resonance information solution. Not surprisingly, there were 
exactly 60 independent significant eigenvectors in the covariance matrix 
of the combined solution minus GEM-T1 corresponding to the number of 
independent (lumped) observations introduced into the solution. But 
only 52 of the 60 eigenvectors were found to be well determined because 
8 represented independent resonance effects determined from very similar 
orbits. For the 52 significant eigenvectors, a satisfactory calibration 
was achieved. The derived eigenvector calibration factor of 1.09 was 
based upon these 52 independent points and agrees well with the value of 
1.14 presented earlier. 
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TABLE B 

GEM-T1 EIGENVALUE PROJECTED ERRORS 
AND CALIBRATION FACTORS 

AVERAGED IN GROUPS OF 50 

CASE A: GEM-T1 vs GEtt-TI without 4-Lasers data 

CASE B: GEM-T1 vs. GEh-T 1 ♦ Surface Gravity ♦ Seasat 

Altimeter 



| GEM 

-T. , 

Kj 

Ki 



< Ah> 

CALIBRATION 

CALIBRATION 

EIGENVALUE 

AVERAGE 

SCALEDAVERAGE 

FACTORS 

FACTORS 

GROUP 

EIGENVALUE 

ERROR(cm) 

CASE A 

CASE B 

1 

58E-20 

0 05 

1.01 

2 12 

2 

66E- 19 

0.16 

0.89 

1.19 

3 

20E-18 

0.29 

0 83 

1.65 

4 

-46E- 1 8 

0.43 

1.14 

1.58 

5 

89E-18 

0 60 

1.09 

1 40 

6 

15E- 1 7 

0 79 

0 85 

1.56 

7 

25E-17 

1.01 

0.93 

1 6^ 

8 

39E-17 

1.25 

1 00 

1 32 

9 

59E- 1 7 

1.55 

116 

1.48 

10 

90E- 17 

1.91 

0 99 

1.16 

1 1 

. 1 3E- 1 6 

2 31 

1.17 

1.25 

12 

17E-16 

2.66 

0 84 

1 .01 

13 

22E-16 

3 01 

1 09 

0 92 

14 

27E- 16 

3.29 

1.14 

0.77 

15 

30E-16 

3.47 

1.12 

0.89 

16 

32E- 16 

3.62 

0.96 

0.83 

17 

35E-16 

3.79 

1.12 

0.94 

18 

39E- 16 

3.97 

0.98 

1.10 

19 

43E- 16 

4 16 

0.78 

1 23 

20 

48E- 16 

4.39 

0.89 

0.93 

21 

.53E-16 

4.63 

1.05 

0 83 

22 

60E- 16 

4 92 

0.93 

0 81 

23 

.69E-16 

5.28 

0.96 

0 91 

24 

82E- 16 

5.75 

1.03 

0.93 

25 

10E- 15 

6.47 

0.88 

0.85 

26 

I3E-15 

7.34 

0.98 

0.83 

27 

.18E-15 

8.47 

0.79 

0.74 
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We have similarly performed these calibrations using the models 
which contained a factor of 10 bias introduced in the data weighting. 
The eigenvector projections yielded a comparable indication of a poorly 
calibrated model (see Table 9). 

There clearly is a good agreement between the two sets of results, 
indicating that both the variances and covariances of the GEM-T1 
solution are well understood. Again, with the exception of LAGEOS’ 
subset solution, all calibration factors shown are within 10# of the 
ideal value of 1 . 

6.4 Summary of Eigenvector Calibrations 

The overall calibration for the eigenvector analysis is expected 
to agree somewhat with that obtained on the basis of the variances of 
the geopotential coefficients themselves. This is because of the 
invariance of the trace (sum of the , variances) of the variance- 
covariance matrix under rotation. This is not surprising because for an 
overall calibration factor obtained from the complete set of coefficient 
differences it holds that: 

T T 

AC AC Ay Ay 

= Tr(AV) = 7 X 

L n 

and it is thus independent of the technique used to determine it. The 
individual calibration factors 


show significant variation among the different data subsets utilized and 
have provided a very interesting spectral decomposition of the errors. 
While it is important to know which parts of the field calibrate best 
and poorest from a statistical point of view, for our purposes it is 
even more important that the parts of the field which the TOPEX orbit is 
most sensitive to are those calibrated satisfactorily. 
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TABLE 9 


1 


SUMMARY OF SOLUTION CALIBRATION FACTORS 
FROM GEM-TI FIELD ASSESSMENTS 


RMS WEIGHTED 
PROJECTED 
EIGENVECTOR 



COEFFICIENT 

CALIBRATION 


CALIBRATION 

FACTOR ONTO 


FACTOR 

GEM-TI 

(GEM-TI) vs. (GEM-TI minus DATA SUBSET) 


4-LASERS (GEOS 1,2,3, BE-C) 

1.06 

0.94 

STARLETTE LASER 

1.10 

0.99 

OSCAR ♦ SEASAT DOPPLER 

1.09 

1.07 

OPTICAL (1 1 SATS) 

0 84 

0 89 

LAGEOS LASER 

1.45 

1.59 

GEM-TI vs. GEM-TI ♦ SURFACE GRAVITY 

0.95 

0 92 

GEM-TI vs GEM-TI ♦ SURFACE GRAVITY ♦ 



SEASAT ALTIMETRY 

0.94 

0 89 

GEM-TI vs SURFACE GRAVITY ♦ SEASAT 



♦ ALTIM 

0.99 

0.90 

GEM-TI minus LAGEOS vs. LAGEOS ♦ 



SURFACE GRAVITY ♦ SEASAT ALTIMETRY 0.95 

0.88 

GEM-T 1 vs. GEM-T 1 ♦ Lumped Resonance 

1.00 

1.06 

Data 



GEM-TI with 10 times the Data Weight vs. 
GEM-TI minus 4-LASERS with 

2 75 

2 45 


10 times the Date Weight 
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Returning to the original question pertaining to the influence of 
collocation on the correlation of the model, our eigenvector analysis 
has shown that there has been no significant aliasing of the coefficient 
variances in the region of the field most important in orbit determina- 
tion. Apparently, collocation controls only weakly determined parameters 
and has little or no effect on those which are well sensed by the data. 
This insight into the behavior of the field could not be gained without 
a detailed calibration of individual coefficients and their corres- 
ponding eigenvectors. 



VII. DIRECT CALIBRATION OF GEM-T1 WITH 1071 5°x5° MEAN 
GRAVITY ANOMALIES FROM ALTIMETRY 


The purpose of this phase of our analysis is to perform a 
calibration of GEM-Tl's errors by comparing gravimetric observations 
directly with corresponding values computed from the GEM-T1 model and 
equating these residuals with the errors produced by the model scaled by 
a calibration factor (k). Previously, except for the direct satellite 
resonant calibration, all other tests performed herein were based on 
subtracting coefficients between two solutions which differed by a given 
subset of data. This difference was then equated to a scaled (k) 
difference of the error estimates (that is, to be more precise, the 
square root of the difference of the error variances). These previous 
calibrations therefore indirectly test the fields against their data. A 
possible concern with this indirect method of data calibration is that 
both geopotential solutions may be affected similarly by unmodeled 
effects and hence, these errors go undetected when forming the field 
difference. The direct calibration of gravity models with independent 
observations as undertaken here will avoid this pitfall. 

Mean gravity anomalies (Ag) for 5°x5° blocks are somewhat 
commensurate in field resolution with the 36x36 harmonic model of 
GEM-T1 . These independent data were used to calibrate the model. 
Although commensurate, analysis of the truncation effects beyond 
degree 36 was considered for the 5°x5° mean anomaly Ag and this effect 
was found to be important. The 5° mean anomalies from altimetry have 
accuracies which are usually better than 1 milligal which is far better 
than the errors from the GEM-T 1 model. The uncertainty ascribed to the 
altimeter gravity anomalies are presented in Figure 12, and are compared 
to the GEM-T 1 model in Figure 13. These figures indicate that the large 
GEM-T1 errors are due to the effects of the higher degree terms which 
are less well known in satellite models. Therefore, the calibration 
will primarily test these harmonic terms. 
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# OF 5 DEG BLOCKS 



7.1 Mean Gravity Anomaly Data for 5°x5° Blocks from Altimetry 

Mean gravity anomalies Ag for 5°x5° blocks were obtained from a 
set of areal mean 1°x1° anomalies, which were derived from SEASAT and 
GEOS-3 altimetry by Rapp (1986). A set of 1071 5°x5° equal-angle Ag 
blocks were formed where each block had a complete set of 25 measured 
1°x1° values. The set covered the entire ocean area and this type of 
data has been used previously to calibrate the earlier GEM models (Lerch 
et al., 1985b). 


7.2 Initial Method of Calibration Excluding Truncation Effects 

Because of the eommensurability of observed 5°x5° anomaly blocks 
(Ag) with those computed from the 36x36 model of GEM-T1 (Ag ) the 
effects of the truncation of higher degree terms were not included in 
the initial field calibrations. Values for the mean anomalies from GEM- 
T1 were obtained by averaging 1°x1° mean anomalies (Ag c ) to form 5°x5° 
mean anomaly blocks. The Ag^ (1°x1°) were computed (Heiskanen and 
Moritz, 1967) from harmonics of GEM-T1 evaluated at the center point of 
the 1°x1° block as follows: 


Ag e ( 1 °x1°) = y 


36 n 

I I 

n=2 m=0 


( £ ) <n-1) 


P m (sin 4>) 
n 


(7.1) 

[ C nm cos mX + S nm sin mX 1 

where r is evaluated on the reference ellipsoid with y = GM/r 2 . The 
5°x5° mean anomaly blocks were formed as: 


Ag c 


25 

l % ( 10x10) 

i = 1 c 1 


25 
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The calibration factor k is related to the mean gravity anomaly 
residual, Ag - Ag c and errors as follows: 


E(Ag - Ag n ) 2 = k 2 a 2 (Ag ) + a 2 (Ag) 


(7.2) 


where a(Ag) and a(Ag c ) are the errors in the mean anomaly data and the 
GEM-T1 errors respectively which are derived below. A global scale 
factor, k, can be obtained from (7.2) by averaging over the 1071 ocean 
points , namely 

M[( Ag - A? ) 2 ] = k 2 M[a 2 (Ag ) ] + M[ 0 2 (aI)] 

C V 

where 


k 


M[Ag - Ag Q ) 2 ] - M[a 2 (Ag)] 
M[ a 2 (Ag Q ) ] 


1/2 


(7.3) 


and M denotes the mean square value. 

The error for each 5°x5° mean anomaly a(Ag) was obtained from 
the 1°x1° errors given with the data as follows: 


o( Ag) 


25 a (Ag 10x1°). 
i= 1 25 


1/2 


/ 5 


This assumes the errors for 1°x1° anomalies are random and 
uncorrelated. These errors are plotted in Figure 13 along with the 
GEM-T1 errors a(Ag c ). The altimeter derived gravity anomaly errors are 
seen to be quite small (RMS = 1 mgal) as compared to the ones derived 
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from GEM-T1 (RMS = 4.5 mgal). For comparison purposes, a second global 
scale factor will be computed from (7.3) where the altimeter anomaly 
data errors are set to zero which will yield a conservative estimate of 
k. This result also shows that these relatively small altimeter block 
data errors are not very significant in the determination of k. 

The calculation of the commission error a(Ag ) from GEM-T1 is 
now described. The mean anomaly errors were obtained from the GEM-T1 
error covariance matrix (V) with the use of a spectral (Pellinen) 
smoothing operator e n (Rapp and Jekeli, 1980). The operator enables 
errors for point values to be averaged over a capsize with area 
corresponding to a 5°x5° block. The error for Ag^ due to coefficient 
errors in GEM-T1 is given as: 

36 n an 

5(Ag ) = Y l l 6 [~) P®(sin 4)) [AC cos mX + AS sin mX] 
n=2 m=0 r n nm nm 


e A AC , 


then 


E E 6 2 ( Ag o ) ] = o 2 (ig c ) = AVA T (7.1)) 

which is evaluated at the center point of the 5°x5° block. For a 
conservative estimate of the calibration factor, the 6 n were computed 
using an area corresponding to the size of an equal-angle block at the 
equator. A more rigorous computation would account for the smaller 

5°x5° areas remote from the equator. 

The importance of the smoothing operator 6 n may be seen by 

comparing o(Ag ) with point values of o(Ag ) where a =1. The 

^ ^ c o 

histogram (Figure 13) for a(Ag c ) shows an rms value of 4.5 mgals, 

which is the square root of the term being scaled by k 2 for the global 

2 — p 

calibration, namely M[o (Ag )]. For point values, M[o (Ag )] 

c c 
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p 

equals about 6.0 mgals . Hence, the effect of B n provides for a more 
realistic calibration factor by the ratio of (6. 0/4. 5) or 1.33. This is 
still somewhat conservative since the smoothing for 6^ is based upon 
equatorial areas (which yields both smaller and a(Ag c )). 


7.3 Calibration Results 

From the global calibration given by (7.3), the value of k 
computed was 

k = 1.1 . 

A histogram (Figure 14) using the global k value was obtained from 
individual calibration factors k^ for each mean anomaly as follows: 


k. 

l 


l A S - A 8 C I i 

7 277-^ 72 2, - , 1 1/2 

a (Ag) + k a Ug e ) | 


(7.5) 


Using a 3 sigma edit of 3.3 (based upon the global value of k = 1.1) 
8 points with large residuals were edited yielding a new global 
calibration factor of: 


k = 1.04 


With use of this k an rms of k ^ from (7.5) gave 


l 4 


1/2 


1062 


1.03 


(7.6) 


(7.7) 
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The value for k rose to 1.07 when the data errors a(Ag) were set to 
zero, showing as predicted, that these errors are not very significant 
in this calibration; the commission error in GEM-T1 is much larger than 
these data errors. 

The geographical distribution of the k^ calibration factors are 
shown in Figure 15 with the 8 points selected for editing shown in 
black. These edited points overlie areas of oceanic trenches. If the 
editing is tightened to 2 sigma, 34 points (Figure 16) are selected for 
editing using (2k = 2.2) as a residual cutoff. These blocks are all 
located in areas of steep geoidal variation, leading us to conclude that 
GEM-TI's truncation error is significant and the effect should be 
included in the calibration. 


7.4 GEM-T1 Calibration with Adjustment for Truncation Effects 

As a result of the relatively large values in the tail of the 
histogram (Figure 14), it was decided to adjust the gravity anomaly data 
for truncation effects, Ag^, for terms beyond degree 36 in the harmonic 
expansion. These new calibrations used altimeter anomalies, Ag ft , 
defined as: 

Ag A (loxio) = Ag( 1 °x1 0 ) - Ag T (10 X 10) (7.8) 

The mean anomalies over 5°x5° blocks, corrected for truncation effects, 
were computed as: 

Ag A = Ag - Ag T (7.9) 

Using the gravitational model of Rapp and Cruz ( 1986a), which is 
complete in spherical harmonics through degree and order 360, a grid of 
\° x Ag T values were computed for terms beyond degree 36. These 
Ag T values were areally averaged to produce Ag T (1°x1°) anomaly 

corrections . 
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Using the adjusted 5°x5° mean anomalies Ag ft in place of Ag in 
(7.9), the calibration process was repeated. The new global value of k 
from (7.3) gave 

k = .87 (7.10) 

which was based upon all 1071 adjusted anomalies with no editing. The 
new k ^ from (7.5) corresponding to each adjusted mean anomaly (Ag ft )^ 
was plotted alongside the previous values in Figure 14. Note the tail 
of the new histogram shows significant improvement, with many fewer 
instances of large disagreements (which required editing for the 
previous calibration). Also, accounting for truncation improved the 
overall calibration factor indicating conservative estimates published 
for the commission error in GEM-T1 . 


7.5 Improved Comparisons of Gravity Models with Adjusted Gravity 
Anomalies 

Based on this new set of 5° anomalies which have now been 
corrected for truncation effects (Ag T ), we have revised some earlier 
comparisons using older GEM models, applying these new anomalies 
instead. For each field the mean square residuals 

l - A 8 0 )2 i 

M(Ag ) = — 

H 1071 

are computed and compared to M( Ag) for the original uncorrected 
anomaly data. Results are presented in Table 10 for comparison showing 
significant reductions in the mean square residuals. 
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7.6 Surface Gravimetry Calibration Conclusions 


An oceanwide set of 1°x1° altimeter-derived gravity anomalies, 
when averaged to form 5°x5° mean gravity anomalies, were found to be 
quite productive in calibrating the GEM-T1 model. The truncation effect 
for 5°x5° mean anomalies was found to be an important contaminant in 
estimating the calibration factor for the 36x36 model. Estimations for 
the global calibration factor reduced k from 1.11 to 0.84 by applying 
estimates of the truncation effect for terms above degree 36 obtained 
from a high degree and order gravity model. This truncation effect, 
when used to correct the 1071 5°x5° mean gravity anomalies, gave a 
significantly improved anomaly data set for testing 36x36 geopotential 
models. The commission error in the field is better tested using these 
new anomaly values, which is the chief concern in field calibrations. 
These adjusted anomalies reduced the mean square residual misclosure 
between altimeter anomalies and GEM-T1 from 29 to 15 mgals 2 . 

Because of the greater sensitivity of the gravity anomaly to the 
higher degree harmonic terms, the calibration results obtained herein 
apply more directly to the higher degree (above degree 10) terms of 
GEM-T1 . 
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Table 10. Improved Field Comparisons of 1071 5°x5° Mean 
Gravity Anomalies Based Upon Mean Square Residuals (MSR) Using 
Original Data (Ag) and the Data Adjusted for 
Truncation Effect (Ag ft ) 


Field 
GEM-T1 

COMBINED MODEL 14 5 

(GEM-T1 + SURF. GRAV. + 

ALTIMETER) 

GEM-L2 40 29 


MSR(Ag) MSR(Ag A ) 

25 mgal 2 15 mgal 2 
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Table 10. Improved Field Comparisons of 1071 5°x5° Mean 
Gravity Anomalies Based Upon Mean Square Residuals (MSR) Using 
Original Data (Ag) and the Data Adjusted for 
Truncation Effect (Ag ft ) 
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MSR(Ag) MSR(Ag A ) 
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VIII. ESTIMATED TOPEX ORBITAL ERRORS FROM GEOPOTENTIAL SOURCES 


Figure 1.2 has indicated that there is a significant improvement 
in the level of coefficient accuracy which has been achieved with the 
development of GEM-T1. It is of interest to assess how these 
improvements translate into radial uncertainty performance on the TOPEX 
nominal orbit. This assessment is now possible since we have demon- 
strated in Section VI that the covariance of the GEM-T1 solution is 
well calibrated and gives a reliable estimate of model uncertainty in 
the presence of model correlation. 

To make these calculations, the GSFC ERODYN Program (Englar et al, 
197S) was used. This program is capable of propagating the full gravity 
model covariance error statistically into an rss position error of the 
satellite's trajectory as a function of time using the variational 
equations involving the partial derivatives of position with respect to 
the force model. The covariance matrices for GEM-L2 (Lerch et al., 
1982), GEM-T1 and a combination solution containing GEM-T1 and surface 
gravimetry /altimetry were individually assessed. In this section we 
give only the radial error contributions from the gravitational field 
alone and consider no other sources of error such as those from the 
tracking systems. 

Three days worth of Doppler data were simulated from a global 
network of 40 stations to provide observations for this analysis. These 
'perfect' observations were made without any consideration for tracking 
system errors and merely reflect the likely data distribution for TOPEX 
from one of the possible tracking scenarios. These observations were 
based on our best knowledge of the nominal TOPEX orbital characteristics 
and were reduced in our standard mode to yield gravitational normal 
equations and the variational matrix of force model derivatives as they 
are computed within our orbit determination program, GEODYN (Martin 
et al., 1987). The covariance matrix of the geopotential solutions 
completed the information required by ERODYN to give a picture of the 
TOPEX orbit errors in time over this tested three day interval arising 
from geopotential uncertainty. 
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TABLE 11 


Radial Orbital Errors (RMS) 
for Three Day Arc Lengths 
Using Calibrated Covariance Matrices 


Geopotential 

Radial 

Model 

RMS Error (cm) 

GEM-L2 

65 

GEM-T1 

25 

GEM-T1 + 


Surface Gravimetry + 
Altimetry 

17 



Table 11 summarizes these results. Shown in this table is the rms 
radial error due only to geopotential effects for the TOPEX orbit over 
this three-day arc length. The results displayed as a function of time, 
are presented in Figure 17. They indicate that GEM-T1 has made 
considerable progress in reducing the likely gravitational modeling 
problems for TOPEX. And with the further addition of alt imeter/sur face 
gravity data types, these errors should be substantially reduced still 
further. While they are preliminary, these findings are grounds for 
cautious optimism. 

Figure 18 shows a breakdown of the GEM-T1 radial errors on TOPEX 
for terms of the same degree and those of the same order. The strongest 
signal is seen from an evaluation of the geopotential error contribu- 
tions by order, where two very significant peaks are found. These high 
error sources are due to the m=1 and m=13 terms. The m=13 harmonics are 
those which have primary resonance with TOPEX and even a small amount of 
TOPEX data will be capable of resolving them to a high level of accuracy 
for they act at a limited and very narrow frequency band. Likewise, the 
m= 1 effects are largely due to the so-called "m-daily ,, perturbations 
which give rise to perturbations with a period of nearly one cycle per 
day for m=1 terms. Again, limited amounts of TOPEX data can 
effectively resolve these "lumped" harmonic effects. Therefore, while 
the overall performance of the fields seems to be improving, this 
preliminary assessment of the most "TOPEX-orbit" sensitive part of the 
error spectrum indicates that limited amounts of TOPEX tracking data (a 
few month's worth) should yield a satisfactory orbit (radial 
errors ~10 cm) even with current models shortly after TOPEX flies. 
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RMS(CM.) GEM-L2 = 66, GEM-T1 = 25, PGS3163 = 19 
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IX. SUMMARY 


This report documents the derivation of calibration techniques 
using independent and dependent data sets incorporating both variance 
and eigenvector analyses. It details their application in assessing the 
accuracy and calibration of the error covariances of the GEM-T1 
gravitational model. Taking full advantage of the "super-computing" 
environment available at NASA/Goddard Space Flight Center, many 
solutions have been compared providing a completeness of field testing 
heretofore impossible within earlier computing environments. The 
results show a model remarkably consistent in stability for the 
calibration of its errors. With the exception of a few known and 
understood high-order resonance terms (and the limitations of the high 
altitude Lageos satellite providing data suitable for the calibration of 
a full 36x36 field), the calibrations show a stability in error 
assessment at the 1 0% level for each of the major data subsets employed 
in this evaluation. The published coefficient uncertainties for GEM-T1 
and its error covariance matrix are herein found to be reasonably well 
calibrated and reliable. For example, the average calibration factor 
(k) for GEM-T1 using nine major sets of data in Table 9 (excluding the 
anomalous result for LAGEOS data) gave k=0.99 (± .08) for the 
coefficient calibration and k=0.95 (± .09) for the eigenvector 
calibration. This is a gratifying result, particularly, since formal 
least squares error formulae based on random variables were employed 
with compensating downweighting factors to account for more general 
formulae involving error sources with unknown systematic effects. 
Appendix A analyzes the mathematical validity of the error estimation 
techique for the gravity model, develops an optimal weighting technique 
with internal self-calibration of the error model, and shows that GEM-T1 
approximates this process. 
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APPENDIX A 

F.J. Lerch and E.C. Pavlis 

WEIGHTING, ERROR COVARIANCE, AND 
OPTIMIZED CALIBRATION FOR BIASED DATA IN GEM-T1 


The least- squares or modified least -squares linearized adjustment 
methods are commonly used to estimate a set of parameters (the solution) 
from a linear combination of observation residuals. The residuals 
(difference between the observations and the computed observables from 
a model) consist of, in addition to the errors arising from the modeled 
parameters, a systematic effect due to unadjusted parameters and a 
random noise component. While the exact cause of these systematic 
effects (biases) is unknown, we will show that they can be accommodated 
for in the solution statistics through the proper weighting of the 
data. It is shown that the ensemble of biases in the GEM-T1 observation 
residuals also have some form of randomness so that the linear 
combination of the data errors tend to average towards zero; this is 
necessary for the signal to be recovered well in the solution. The 
analysis in this appendix is an attempt to (1) characterize the bias 
errors, (2) weight their effects along with the noise errors, 
(3) develop an approximate form of the error covariance, and (4) provide 
an optimized technique for weighting the various data subsets. The 
calibration of the error estimates for our solutions hinges on deriving 
the optimal data weights; we use a combination of solutions which use 
data subsets in these assessments. 

The estimation technique for the GEM-T1 solution is a weighted 
least squares process with a priori constraints on the size of the 
adjusting coefficients. The inverse of the weighted least -squares 
normal equations (N) is used as the error covariance of the estimated 
parameters. This process is common in least-squares solutions in a 
formal sense, and it is formally correct when the mathematical model is 
complete and thus the computed residuals are unbiased. Since the GEM-T1 
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observation residuals display a biased behavior due to unmodeled 

effects, a special weighting of the normal equations is employed to make 

1 

the simplified error estimation (N ) applicable as an approximation for 
the error estimate of the solution. An optimization of this weighting 
process is developed herein which automatically provides for a well 
calibrated error estimate of the gravity parameters. The appendix is 
presented in three parts: 

A1 . Weighting for Biased Data 

A2. Error Covariance for Biased Data Reduced to Simplified Form 
( N - 1 ) , and 

A3. Optimization of Weighting With Error Calibration 
There is also a Summary section to the appendices. 


A1 . WEIGHTING FOR BIASED DATA 
A1 . 1 Experimental Analysis 

The characteristics of the weighting system are based upon the 
actual performance of the observation residuals. The tracking data is 
post-processed, in this case using GEM-T1, over orbital arc lengths used 
in the GEM-T1 solution, and the data residuals are analyzed to determine 
the systematic characteristics seen in the passes of tracking data. Let 
the residuals r be denoted as 

r = r g + n (A 1-1) 

where r s is the systematic part of the residual and n is the random 
noise. Typically the systematic effects appear as a straight line (when 
the residual is plotted against the observation derivative in time). For 
example, for laser ranges (R) over a pass of NOBS points 
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( A1-2) 


r = r + R x 
s o 


where r Q is an offset constant, t is a timing constant, and R is the 
range rate. Let the average data residual per pass represent a bias 


Er 

Er _ s_ En 

b = NOBS " NOBS + NOBS 


(A1-3) 


where the symbol (e) indicates summation over the components in the 
vectors (r, r , and n). 


Since the random noise will largely average out over the typical number 
of points seen in a pass, the significant effect is given by 


E(b) 


Er 


NOBS 


i 0. 


(A1-4) 


Note that for symmetric coverage about R = 0, (point of closest 
approach) : 

E(b) - r Q . 

In the analysis presented here using a sample of passes (Brown, 
1988) fit with GEM-T1 , it was found that the biases (as in (A1-3)) are 
the dominant part of the systematic effects seen in r s (from ( A 1 -2 ) ) 
which may be associated with a bias in the measurements. For example, 
on 450 passes of laser tracking acquired on GE0S-1 during 1977 and 1978 
(mostly SAO/NASA tracking), the rms of the residuals, o(r), and the 
biases a(b) for all the passes were 

o(r) = = 70 cm 

a(r ) = 60 cm 
s 

a(b) = 50 cm 

NOBS = 60 average pts. per pass (after data sampling) 
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where o, refers to the residuals for a specific satellite data type t 

L/ 

(i.e., GEOS-1 laser ranges). 

The histogram of the biases for the 450 passes appear as a normal 
distribution with an average bias 

b = 0.2 cm. 

Inspection of the biases (per pass) in time plots on a daily basis show 
fluctuations that appear random. At this point it is unimportant to 
characterize whether these biases are in the data themselves or whether 
they are a result of unmodeled (or poorly modeled) force modeling 
effects. What is important is that these errors yield data residuals 
which are quite systematic over a pass and are describable using the 
model given in (A1-2). 

During the 1 980 period the laser data on GEOS-1 was not used in 
GEM-T1 , but a similar analysis of these data gave an improved set of 
residuals: a = 35 and o(b) = 25. This result is typical for the laser 
systems, especially for the SAO systems which were improved over this 
time interval. Similar analysis on a pass -by -pass basis is being 
developed for other satellites. Preliminary results on Starlette and 
Lageos indicate that much smaller and less biased residuals are 
obtained. Because of the significant improvement in the residuals and 
biases on GEOS-1 laser data between 1977 vs. 1980 with the same gravity 
model, it would appear that a significant source of the biases for the 
earlier time frame is the measurements themselves (or related pre- 
processing) . 
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A1 .2 Effects of Biased Data on Least Squares Estimation 


To examine the effect of biased data (unmodeled effects) on the 
solution, a single pass of tracking data with residuals r contributes to 
the overall weighted least-squares solution (as in GEM-T1 ) , as 

Q = f r Wr N points for the pass ( A 1 -5 ) 

where 

W = h/°Z] diagonal weight matrix 

a = rms of residuals for satellite data type t 

f = scale factor to be estimated to account for effects of 

biased data in the solution 

Let the bias* be constant per pass and random from pass to pass with 
variance for a given data type as 

E(b 2 ) = o 2 (b) 

Denote the residual r as 

r = Ax + b + n 

= r + n (A1-6) 


where 


x = gravity and other unknowns 


*The bias b and other constants are used interchangeably as a constant 
and a vector constant. 
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n = random noise 


b = bias 

r = systematic effect as in (A1-2). 
s 

For simplicity, consider the systematic part of r from (A1-6) as 
r^ = Ax + b 

o 

Ax = g + Ag (A1-7) 


where 


zA. x 
- 1 

g = NOBS 


average offset (i = 1 to NOBS) 


( A1-8) 


and let g be estimated from the biased data. Minimizing Q wrt g will 
give the normal equations as 

f G T WG g = f G T Wr (A1-9) 

M = f G^WG normal matrix (A 1-10) 

where G T = (1,1,1...1) NOBS ones for the pass. From (A1-6 through A 1-9) 


- _ zr ZAx z(n+b) 
g ‘ NOBS = NOBS + NOBS 

- , zn 

- g + b + N0BS 

with error and variance respectively as 


g 


- g = b + 


zn 

NOBS 
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(A1-11) 


o 2 (g) 


a 2 (b) + 


o 2 (n) 

NOBS 


From ( A 1 — 1 0 ) the least-squares error estimate of the variance is 


o 2 (g) = M 



f NOBS 


(A1-12) 


Equating (A1-11) and (A1-12) gives the scale factor 


f = 


a 2 (n) 


2 



+ NOBS a 2 (b) 




(A1-13) 


( A 1 — 1 4 ) 


2 

since the effect of the noise is small compared to NOBSxa (b). 

p 

Noting that a (b) varies with the data type 


a fc (b) = a(b) 


and since it is a significant part of the rms of the residuals, 
then 


a(b) 


a (r) 

— ~ 1 

a t (b) 


(A1-15) 
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Hence (A1-14) becomes, using (A1-15) 


~ NOBS 


(A1-16) 


which is a factor used to downweight the data to account for the 
unmodeled bias in the pass. 


For GEM-T1, f=.02 was found to be optimal; this corresponds to N=5Q on 
average which agrees well with the overall distribution of GEM-TI's 
data. This is especially true for the sampling used for the SAO laser 
data where a typical pass contained 50 to 60 observations (Brown, 1988). 


A 1 . 3 Optimal Estimation By Adjustment of Bias Per Pass of Data 

Herein, we will analyze the effect of the bias adjustment on the 
weighting and error estimation from a pass of data as defined in Section 
A1.2. From (A1-6) the residuals r are unadjusted for bias. Designate 
the adjusted residuals for bias as 

r=r-b=n+Ax+b-b 
where as before 

r = n + b +Ax 

An a priori constraint a ^ = a(b) is applied to the bias in the 
adjustment. Using the symbols and definitions of Section A 1 . 2 we first 
minimize (similar to equation ( A 1 —5 ) ) wrt b, solve for an a priori 
value of b, and transform by back substitution for b as follows: 

n ~T „ "2 . 2 

Q, = r Wr + b/o A 
b A 
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1 T A b 

1 Q b = -G 1 W r ♦ T = o 

°A 


b = + 


z r 


NOBS+q 


q = (a fc /o A )‘ 


_ r T VJ r (N+q-1 ) z r 2 - (z r) 2 
b NOBS+q 2 (N0BS+q) 

U 


( A 1 — 1 7 ) 


Note that if the residuals are dominated by the bias b, then 
o t (r) - or q=1. In this case, with q=1 and 


r = constant bias over the pass, then 


o - mix 

w b NOBS+1 ‘ 


( A 1 — 1 8 ) 


This result is very similar to the weighting in ( A 1 -5 ) where by ( A 1 — 1 5 ) 
and (A 1-16) 

Q - Q b 

with 

f “ NOBS 

For GEM-T1, f=.02 corresponding to NOBS=50 on average. Hence the 
weighting in GEM-T1 is near optimum for the case where the laser 
residuals are dominated by a bias and all passes would ideally have 50 
points . 


If, in addition to a bias, we consider a timing error t (as in 
(A 1-2)) in the residuals with balanced sampling about R=0 (i.e., a pass 
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which is symmetric about the point of closest approach of the satellite 
to the station) where 

r = r = b + R t 
s 

then (with q=1) in (A1-17) 

r T W r z(R t) 2 ( A 1 - 1 9 ) 

y b “ NOBS+1 + 2 

°t 

This case shows that the downweighting factor f only principally 
influences the bias. The approach would still be optimum if the timing 
errors were largely attributed to the gravity field signal. 

Therefore, this method is ideal when biases dominate the 

measurement residuals. It has great benefit in the proper relative 

weighting of the data, particularly when there is significant 

variability in the number of points per pass and in the magnitude of the 

biases a (b). The bias variability can be systematic chronologically as 
t 

will happen as the laser systems are improved. Of most benefit, this 
approach reveals the reduction of the adjusted residuals o fc (r) when 
biases are removed as compared to the a fc (r) when biases are not 
removed. Based upon this approach for error estimation, a reduction of 
the ratio of o fc (r)/o t (r) should occur in the errors (which is a 
significant reduction for laser tracking when normal points are 
employed) since o t (r) would be small. 

It should be noted that A1.2 accounts for the effect of bias on 
the solution when it is not modeled whereas A1 .3 considers its effect 
when it is modeled, which is the recommended approach. Since our 
solution (GEM-T1 ) does not adjust for biases (per pass), this situation 
is analyzed in the subsequent sections. 
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A2. ERROR COVARIANCE FOR BIASED DATA REDUCED TO SIMPLIFIED FORM (N -1 ) 
A2.1 FORMAL LEAST -SQUARES ERROR ESTIMATE 

The weighted least-squares solution of the normal equations with 
random observation errors (n) and complete modeling of parameters in the 
solution (X) gives a simplified form of the error covariance matrix 
V(x). A more general form is given when the observation or modeling 
errors (e) have systematic effects (e g ) in addition to the pure random 
errors attributable to (n). 

The following development is given for the error estimation, 
namely : 


r e 0 - C 
e e 0 - C 


= n + e 


s 


X 



C = C + A x 


observed minus computed quantities from solution X 

C computed from true values X 

n - random noise 
e a - systematic error 

error in solution (A2-1) 

A is matrix of partials for Taylor's 
expansion (linear) about the true X. 


From the above 

A 

r = 0- C- Ax 


= e - A x 


( A2-1a) 


The least-squares weighted normal equations with weight matrix W (to be 
defined) become from (A2-1, la). 
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( A2-2) 


A T Wr = A T We - A T W A x = G 
with error 


x = (A T W A)" 1 A T We ( A2-3) 

Note that if e is biased, then x is expected to be biased, as is the 
case within our solution. The error matrix for the solution is 

V(x) = E(x x T ) = E[x-E(x) ] [x-E(x) + E(x) E(x) 

= (A T WA) _1 A T W E(ee T ) WA^WA)" 1 (A2-4) 

We will refer to V(x) as the error covariance. For the most generalized 
case where e has systematic effects let 

W** 1 = E(ee T ) which is non -diagonal, 


then 


V(x) = (A T WA) _1 = N” 1 ( A2-5) 

This case is not generally useful since W is non-diagonal and generally 
unknown given our lack of knowledge of all the unmodeled errors. For 
the most simplified case let 

e = n 


then 




which is a diagonal matrix 


and 


V(x) = [A T WA] _1 = N 


( A2-6) 
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which is the formal least squares error covariance estimate. 

Our solution for the error x is a mixture of (A2-5) and (A2-6). 
We cannot apply W 1 of (A2-5) to the solution (A2-3) since the weight 
matrix is non-diagonal and unknown. Further we cannot apply (A2-6) 
since the error covariance is assumed to be diagonal and our errors by 
definition include biases which are strongly correlated within a pass of 
tracking data. Hence our errors follow (A2-4) where ideally we model 

W = f [l/o 2 (e>] which, again, is a diagonal matrix. (A2-7) 

Here f is a downweighting factor accounting for the correlated 
systematic error effects in modeling the data anS o^e) is estimated 
from the rms of the residuals of a given data type for a given satellite 
based upon its post-solution performance. It remains for us to show 
that (A2-4) with (A2-7) can be approximated by (A2-6) under empirical 
conditions, albeit, which are reasonably simplified. 


A2.2 EMPIRICAL MODELING OF SYSTEMATIC EFFECTS 


For our case consider systematic effects (e s ) on a P ass by P ass 
basis for a given satellite's data. We typically find these systematic 
signatures to exhibit a straight-line character as follows (see Section 

A 1 . 1 ) : 

(A2-8) 

e = n + e s 


Ae 


= e a (t) - % (t m ] 


where 


e 


s 


straight-line fit 
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b - constant offset centered at the midpoint (t m ) of 

the straight line 

t - constant timing error 

e - time rate of e 0 

s s 

n - random uncorrelated noise 

For each data type the biases over all passes can be characterized as 
(see Section A 1 . 1 ) 

E(b 2 ) = a 2 (b) (A2-10) 

and our experience has shown that it is reasonable to assume these 
systematic trends are uncorrelated among passes (certainly, for 
weighting purposes), then 

E(b. bj) = 0, i * j . (A2-11) 

We will further assume that within a pass the bias b is orthogonal 
to the timing error term 

bZ Ae x = 0 (A2-12) 

s 

which is reasonable (in our case) since the orbits are nearly circular 
and the coverage is usually somewhat symmetric about the midpoint t m of 
the straight line. 
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A2.3 ERROR COVARIANCE AND SIMPLIFIED MODELING 

Let the contributions to the least- squares normal matrices be 
defined as follows: 


(a) for each tracking pass p of observations of a given data type, t, 
with the number of observations equal to NOBS^ 


N . = (A WA) 
p , t P » t 


(A2-13) 


R , = (A A We) n , 
P»t P»t 


where from (A2-7) with f=fp^ t 


W = W . = [f . / a (e) 
p,t 1 p,t t 


which is a diagonal 
matrix of dimension 
NOBS,. 


(b) for each data type t 


N, = E N^ 

t P,t 


(A2-14) 


R, = I R . 
t pt 


(c) for all data types 


R = I R fc (A2-15) 

The normal equations with solution x are 


N x = R 
x = N" 1 R 


(A2-16) 
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and the error covariance is ( a priori constraints to be added later) 


V(x) = E(x x ) 


-1 T -1 

= N E(R R ) N 1 


(A2-17) 


We will show that 


E(R R T ) N" 1 * I 


(A2-18) 


and hence (A2-17) is simplified to 


V(x) = N” 1 


( A2- 18a) 


From (A2-11), (A2-12), (A2-8,-9), and the condition of 

uncorrelated random noise (n) we have 


E I R p,t R q,t’ 1 = 0 or tn ’ 


(A2-19) 


= E f R p,t R p,J P=q and t=t ’ 


hence, using (A2-13) through (A2-16) 


E [R R T ] = E [(Z R ) (Z R ) T ] = E [I R R^ 


( A2-20) 


= E [V y (R R T )] , 
t p P,t P,t 

E [R fc Rj] = E [I (R pit R ; jt )] 

= E[I (A T We) (fl T Wef ] t , 


( A2-21 ) 


and 


(R . R T J = (A T W e e T W A ). 
p,t p,t p p p p p p't 


( A2-22 ) 
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For each pass of data, from (A2-1, -8, -9) 


(0-C) = e = n + e 

p p p s,p 


(0-C) = r = e - A x 

P P P P 


for NOBSp data points. 


(A2-23) 


Since gravity and other parameters modeled in the solution (A x) produce 

~ P 

misclosure errors in the residuals r within a pass of data which (like 

P 

e s ) largely fit to a straight line (see equations A 1 - 1 , -7 , -8 ) , we will 
average the signature A p over the pass (A1-8) for each component k of 
the solution. Thus for each data type let 


A = A 
P 1 P 




MOBS 



(A2-24) 


where 


G T = ( 1 , 1 , 1 ... 1 ) NOBS p ones 

k = 1 to K solution parameters. 

Note that in (A2-12), with b replaced by A^, the orthogonality holds. 

Using (A2-24) and (A2-8, -9, -1 1 ,-12,-21 ) in (A2-22) with A = A , 

P P 

then 
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(A2-25) 


E ( R p,t R p,J =E [A P W P 6 P e pV A p ] t 


= K «P A p't 


where 




2 t \ 

(e) 

+ N0BS P ,t °? (b> ' 


Notice that (A2-26) may be approximated by 



1 


NOBS 

p,t 


(A2-26) 


(A2-27) 


if we assume (as found in Section A 1.1) that the variance of the 
bias, o^(b), dominates the errors o^(e) for a typical pass. (In the 
presence of biases of weighting this type tends to equalize the data by 
passes instead of by number of observations.) 

Using the result (A2-25) in (A2-20) and (A2-21), then with 

A p = A P 

E [R R T ] = II (Ap w A) t (A2-28) 

t P ^ H 

This simply gives from (A2-13, -14, -15) 

E [R R T ] = N (A2-29) 

yielding the desired result (A2-18) and hence 

V(x)=r 1 . (A2-30) 
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This represents the simplified form of the weighted least-squares 
inverse . 

For practical application, the weighting W p>t in (A2-26) can be 
optimally achieved by solving for a bias per pass as indicated in 
Section A1.3 (see equation A 1 — 1 9 ) . The next level of simplification is 

shown in Section A3 of the Appendix. Here an average factor of f p>t in 
(A2-26) is employed which may be represented by 

o 2 Ae) 

f = f = i (A2-31) 

P,t o 2 ( n) + NOBS t c*(b) 

where a mean number of points per pass (N0BS p ) for each data type t is 
employed. In Section A3 the weighting factor f t in actual practice is 
resolved for each data type t from analysis of subset solutions 
associated with each data type requiring the condition (A2-30). 

A2 . 4 MODELING THE A PRIORI KAULA CONSTRAINT 

Before showing further simplification of the weighting used in 
GEM-T1 , it is desirable to introduce the Kaula constraint equations 
which we may regard as data type t=0. The equations for the errors in 
the solution are obtained directly from the constraints on the 
coefficients (C) as follows: 

firstly, denote the coefficients of degree a as 


am 


am 


m = 0 to a (s^ 0 =o), 


the entire set by 


(A2-32) 
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c = [ c J 


the solution errors AC (a subset of x) by 


AC = C - C , 

A 

then the error e Q and residual r Q (data type t=o) are 
e Q = O(zero) - C 

A A 

r = O(zero) - C = e + AC 
o o 

The contributions to the normal equations for t=o become 


N AC = R 
o o 



R o = w o e o 


of (e 0 ) = \ (10- 5 /, 2 ) 2 


(Kaula’s rule). 


(A2-33) 


(A2-34) 


The weighting factor f * 1 can be contrasted to the actual tracking 
data types, where f * .02 was found to be appropriate. From equation 
(A2-27) each coefficient behaves like a bias with one observation point 
as in (A2-33). 

We statistically treat the coefficients C in the error equation 
e Q , (A2-33), similar to the biases for the tracking data, namely 


E (C. C ) = 0 i * j 

E <C 2 ) = o 2 (C.) (A2-35) 
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With this treatment the normal equations (A2-34) with the error e Q (A2- 
33) may be combined as data type t=o in (A2-15, -16, -17) resulting in 
( A2-I8a) . 

A2.5 RANDOMNESS TESTS FOR THE ENSEMBLE OF BIAS PARAMETERS AND 
COEFFICIENTS 

The degree of randomness of the biases and the coefficients (as 
with random noise) is important for the error solution x in (A2-16), 
which is a linear combination of all these error effects. The effect of 
randomness, with proper weighting, permits (x) to average toward zero. 
A measure of the randomness of the biases and the coefficients is given. 


A measure of the randomness of the ensemble of the biases for a 
given data type t is 

IbJ 


Z = 


° (b t ) 


< 2 


(A2-36) 


where 


z b. 
i 


t ' NB. 


- a t (b) 

a(b, ) = Y/2 ’ 

t (NB t ) 


with NB t being the number of passes. 

Similarly, for the randomness of the coefficients 


C 


a 


l (C *m + S vJ 

m 

2 % + 1 


(A2-37) 
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(Modified Kaula's Rule) = o(c^) = 10 *V/2/a 2 
a(C ) = (10" 5 //2/Jl 2 )/(2Jl+1) 1/2 

X/ 



a(C t ) 

Some results are available. 

From Section A1.1, the biases for the GE0S-1 (1977-78) laser data 
employed in GEM-T1 gave (for 450 passes) the following values: 

b, = .2 cm 

X/ 

a (b) = 50 cm 
a(b^.) = 2 cm 
Z = 0.1 

A time ordered plot of the biases is given in Figure A2-1 and a 
randomness test per tracking site for these biases is given in Table 
A2-1 showing strong indication of stochastic behaviour from pass to 
pass . 


Using a solution (PGS-3325) which is a combination of GEM-T1 
normals, surface gravity, and altimeter data, the statistics in (A2-37) 
for the randomness of the coefficients are listed in Table A2-2. Since 
the constraint dominates the higher-degree terms in GEM-T1 , it is 
important that the Z(&) improve with degree a, which is supported by the 
tabulated results. The random behaviour of the coefficients causes a 
minimal error in the solution when the "apparent" biases (C) in (A2-33) 
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Figure A2-1 . Pass Biases 
GEOS-1 Laser Tracking Stations 
(187 passes - 8 months, 1977) 
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TABLE A2-1 

BIAS RANDOMNESS TEST PER LASER SITE 
(Z < 2.5) 


LASER 

STATION 

PASSES 

o(b) 

CENTIMETERS 
b a(b) 

Z 

AREQUI91 

68 

54 

9 

6.5 

1.4 

HOPLAS91 

43 

40 

11 

6.1 

1.8 

ORRORL91 

46 

41 

-15 

6.0 

-2.5 

RAMLAS1 

4 

57 

3 

28.5 

. 1 

STALAS1 

6 

17 

-2 

7.0 

-.3 

NATAL9 1 

15 

98 

-15 

25.3 

-.6 

G0RF1001 

4 

42 

-34 

21.0 

-1.6 

G0RF0652 

1 

177 

177 

— 

— 

AVG 

(187) 

54 

1 

3.9 

.3 
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TABLE A2-2 


Randomness Test for Potential Coefficients* 
of Degree l (Z < 2.5) 

(Coeff. units of 10"^) 


Degree 

6 

No. of 
Coeff. 

RMS 

(Coeff.) 

Kaula 

x 0,7 ....... 

c i 

° ( ? 

Z 

3 

7 

1122.7 

785.7 

807.90 

296.96 

2.72 

4 

9 

529.2 

441.9 

161.19 

147.31 

1.09 

5 

11 

351.4 

282.8 

-105.43 

82.28 

-1.24 

6 

13 

251.3 

196.4 

-158.04 

54.48 

-2.90 

7 

15 

193.2 

144.3 

23.43 

37.26 

0.63 

8 

17 

118.9 

110.5 

25.65 

26.80 

0.96 

9 

19 

97.7 

87.3 

10.50 

20.03 

0.52 

10 

21 

77.7 

70.7 

-26.82 

15.43 

-1.74 

11 

23 

53.5 

58.4 

-19.00 

12.19 

-1.56 

12 

25 

31.2 

49.1 

6.17 

9.82 

0.63 

13 

27 

45.7 

41.8 

8.49 

8.05 

1.05 

14 

29 

26.0 

36.1 

0.89 

6.70 

0.13 

15 

31 

25.3 

31.4 

2.38 

5.64 

0.42 

16 

33 

23.3 

27.6 

-3.87 

4.81 

-0.81 

17 

35 

20.0 

24.5 

-0.06 

4.14 

-0.01 

18 

37 

18.6 

21.8 

-0.84 

3.59 

-0.23 

19 

39 

15.2 

19.6 

0.35 

3.14 

0.11 

20 

41 

13.7 

17.7 

1.06 

2.76 

0.39 

21 

43 

15.1 

16.0 

2.03 

2.45 

0.83 

22 

45 

14.0 

14.6 

1.39 

2.18 

0.64 

23 

47 

11.7 

13.4 

-0.79 

1.95 

-0.41 

24 

49 

10.1 

12.3 

-0.15 

1.75 

-0.09 

25 

51 

11.5 

11.3 

-0.12 

1.58 

-0.08 

26 

53 

8.8 

10.5 

1.23 

1.44 

0.86 

27 

55 

7.8 

9.7 

0.59 

1.31 

0.45 

28 

57 

8.8 

9.0 

0.47 

1.19 

0.40 

29 

59 

7.6 

8.4 

-1.59 

1.09 

-1.45 

30 

61 

8.2 

7.9 

-0.66 

1.01 

-0.65 

31 

63 

8.0 

7.4 

-0.10 

0.93 

-0.11 

32 

65 

6.9 

6.9 

0.06 

0.86 

0.07 

33 

67 

6.9 

6.5 

-0.51 

0.79 

-0.64 

34 

69 

7.2 

6.1 

-0.51 

0.74 

-0.70 

35 

71 

6.7 

5.8 

-0.72 

0.69 

-1.05 

26 

11 

12 

16 _ 

-0.18 

0.64 

-0.28 


* PGS-3326: 50 x 50 model GEM-T1 without Kaula constraint with surface gravimetry and 
SEASAT altimetry. 
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are properly weighted. Hence, any downweighting of f Q from unity in 
(A2-34) or misuse of a^(e Q ) by ig nor i- n S the scale factor of 1/2 will 
allow the errors in the coefficients to increase and falsely represent 
their power. Note that when the Kaula constraint is present, the power 
in the coefficients should become zero for high degree a at the point 
when the satellite tracking data just begins to have negligible effect 
on the solution. 

A2 .6 APPLICATION TO GEM-T 1 

In the GEM-T 1 solution, a more simplified case of the weighting of 

W i_ and f_ f in (A2-26) was employed. The more general weighting in 
P » b P » “ . 

(A2-26) is based upon NOBS p t (points per pass) while the simplified 

weighting is based upon an average NOBS over all passes and data types 

in the solution, namely 

W t = f/o t (e) (A2-38) 

f = 1/NOBS ( A2-39) 

= 1/50 (for GEM-T 1) 

In GEM-T1 , a 2 (e) was approximated basically by the rms of the 

V 

observation residuals from iterative test solutions. It was allowed to 
vary from the rms by use of tests with other data, in order to account 
for the variability in NOBS for different data types (N0BS t ) and other 
variable effects in (A2-26). Furthermore, in the special case of the 
optical data where considerably fewer points per pass were taken than 
the laser data, an effective f = 1/10 was employed by applying a fc (e) to 
1/5 of its corresponding rms residual value so these optical data retain 
some influence on the solution. 

Since the weighting represents an average for the different 
satellite data types and since the simplified form of V(x) in (A2-28a) 
is an approximation, an optimal technique for iterating on a weighting 
factor per data type 
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(A2-40) 


w t = f t / 0 t (e) 

and a scale factor (k) for the error covariance V(x) has been developed 
and is presented in the next section. GEM-T1 followed this procedure in 
somewhat of a piecewise fashion and hence the model may be further 
refined with use of this more rigorous optimal weighting approach. 

The data type t should be considered broadly in general 
applications of this approach. There is variability in the precision 
and accuracy among the different laser tracking sites, hence the data 
type in (A2-40) should reflect this situation. The general improvement 
seen in modern tracking performance as time progresses is another 
consideration. Although the former condition has not been generally 
applied to develop station specific weights in (A2-40) in GEM-T1, the 
latter condition is considered in the individual weighting factors 
applied to the normal equations for orbital arcs of tracking data which 
sequentially span periods of time. In GEM-T2, which will employ some 30 
different satellites, we utilize variable weighting factors for 
satellite tracking data of a similar type. In particular, the laser 
data from GEOS-1 and GE0S-3 spanning 1975 through 1978 is given smaller 
weight than the newer data in 1980 (see Section A1 . 3 ) - 

A2.7 SUMMARY 

If a more generalized description of the different data types is 
utilized for the application of the weighting in (A2-26), the better the 
approximation of the simplified inverse N" in (A2-30) will become. In 
summary, is the inverse of the least-squares weighted normal 
equations and represents, as an approximation, the complete error 
covariance expression of the solution parameters in (A2-17), namely 
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(A2-41) 



V(x) 


Y (A T W A ), 

_tt P p p p 




The result is considered an approximation because it is based upon an 
average signature A^ = A^ over each pass p of tracking data for the 
gravity parameters and other modeled parameters. Furthermore, the 
errors ( e p) » in addition to random noise, were assumed to have 
systematic effects which were empirically modeled as a bias and a timing 
error for each pass of tracking data. The result has also considered 
the effect of the Kaula constraints by modeling it as data type t=0 as 
given in Section A2.4. 
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A3. OPTIMIZATION OF WEIGHTING AND AUTOMATIC CALIBRATION 

There are two major concerns in the proper weighting of least- 
squares normal equations: 

(a) weighting of individual observations corresponding to the 
expected accuracy of the observation residuals (a fc ) for a 
given data type t (c^ is computed from the iterated, rms of 
the observation residuals), and 

(b) the effects of unmodeled biases and forces on the solution by 
applying a downweighting factor (w fc ) for the normal equations 
(see Table 1 for GEM-T1 weights). 


The normal least-squares approach, given c^, accounts for (a) but not 
for the effect of (b) without some special analysis and downweighting 
process. In our approach we obtain on a few test arcs for each data 
type based upon trial (iterative) solutions of the weighted normal 
equations (w ) . It is the purpose of this section to derive the weights 
w fc . We will show that these weights can be obtained in a process that 
automatically calibrates the error estimates of the gravity model. 

The optimized process of determining the weights (w^.) for each 
data type in the solution (with an automatic calibration of the error 
model) is a refinement in the determination of the GEM-T1 solution. As 
before, this process is based upon forming subset solutions and 
comparing these solutions with the complete solution. The subset 
solution is formed by deleting a given satellite data type from the 
complete solution. The weighting and calibration is then established by 
requiring that the difference between the solutions and the error 
estimates be compatible. It will be shown that in the process of 
converging the weights (w^.), one automatically obtains a calibration of 
the error estimates. Hence this process, when applied to GEM-T1 data 
sets, can then be tested by calibrating the new (refined) GEM-T1 
solution with independent surface gravity data by the method given in 
Section VII of this report. 
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Solutions are formed after summing the weighted least-squares 
normal equations. For simplicity let each subset of normal equations 
(N t ) correspond to the reduced normals, where only the gravity 
parameters x remain, and assume all other parameters have been accounted 
for through the back substitution process. Denote the subset normals 
for a given data type by 

A 

w N x = w R t = 0 to T representing all of the data 

L in the solution, 

where t=0 corresponds to the Kaula constraint equations 

p 

with w o =1/o^ (see Section A2.4). 

It is convenient to write 


where a fc is given by (a) above and f t accounts for the downweighting 
effect in (b). Note for the case where biases are the dominating errors 
in the data, equations ( A 1 — 1 5 ) and ( A 1 — 1 6 ) show that f^. is expected to 
correspond somewhat to the reciprocal of the average number of points in 
a pass for the given data type. For the case of GEM-T1 (Section A2) the 
overall value of f=.02 would correspond to an average of 50 points per 
pass for all data types. 

Denoting x as the solution (S) for all data subsets and x t as the 
subset solution (S^.) of all data except data type t, the difference in 
the solutions corresponds to the difference in the errors between the 
solutions, namely 


x-x^ = [X-X(true)] - [X^-X(true)] = X-Xj. 


Hence, for the analysis we will assume as in (A2-1) 
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X = X - X (true) 


x t = X t - X (true) 


The following quantities are defined from the weighted normal 
associated with the solutions S and S fc : 


w t N t x = w t R t 


N t = A t A t 


R t = A t e t 


subset normals for 
data type t 


N = T w . N . 

jft J J 

R = l w , R, 

J*t J J 

N x t = R 

(N + w fc N fc ) x = R + w t R t 

v(x t > = r 1 

v ( x) s + 

x = V(x) (R + w fc R fc ) 

x. = V( x ) R 
t t 

IT. 

V(x fc -x) = E[x fc -x) (x fc -x) ] 


data type t removed 


subset solution 
complete solution S 
error variance of x^ 
error variance of x 


error variance x-x^. 


= V(x fc ) - V(x) 
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equations 

(A3-1) 

(A3-2) 

(A3— 3) 
(A3-4) 


(A3-5) 



(A3-6) 


V(x-x ) = N" 1 - (N + w fc N t ) _1 

* w fc N _1 N fc N" 1 t*0 

Letting Tr denote the trace of a matrix 

Q t = (x t -x) T (x t ~x) = Tr[(x t -x)(x t -x) T ] 

then a scale factor, k t , is defined as 

Q t = k fc E(Q t ) = k fc Tr e[ ( x fc -x) (x t -x) T ] . 

Since k t scales the variances of the errors, V(x-x t ), i 
inversely proportional for scaling the weights. Hence, an 
weight w' is 



This result can be shown more directly from use of (A2-25, 
where 


w t R t = l R P ,t 

p 


= A t w t e t 


w t = V 


then, with use of an updated weight , 


E < R X> = E t A Vt A1 = M t > 

for which, as in (A3-1) 


(A3-7) 


(A3-8) 

will be 
adjusted 


26, -31) 


(A3-9) 
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In order to establish the result we need (x-x fc ) to the same 
approximation as ( A3— 6 ) , namely 

x - x t " w t r 1 R t 

which gives from (A3-5), ( A3— 6 ) , and (A3-9) 


T 

(x-x fc )(x-x fc ) 




N fc N 



V(x-x t ) 


and then from (A3-8) 


(A3- 10) 


k 


t 



(A3-11) 


In the above development all the formulae are based upon the a priori 
weighting (w t ) with the exception of (A3-9) where the updated weight 
(w") is based upon the observation errors e^. including unnmodelled 
effects. Thus, the adjustment of the weights is obtained from x-x t 

which is a linear combination of these errors e^.. 

In summary, the adjusted weight for each data type t can be 
computed from 

V = k 1 w t (A3 ' 12) 

where from (A3-8) 
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(A3— 13) 


1 - 

K t ' Tr V(x -x) 

This process should be iterated until the weights converge for all t 
(k t =1). An overall calibration factor k from (A3- 12) when summed over 
all t is given as 


z Q f 

Lr - - 1 

I Tr V(x.-x) " 
t 

when the weights converge, since k^.= 1 for all t. 


(A3- 14) 
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SUMMARY 


A summary is given for the three areas of analysis in the Appendix 
( A 1 , A2 , and A3). Appendix A1 has shown that the solution scaling factor 
for overall weighting of the data, f=0.02, found to be optimal for 
GEM-T1 , is attributable to the non-random characteristics of the 
tracking data residuals in post-fit analyses of these data. While the 
specific cause of this anomalous data behaviour is not presently known, 
the net result is that each pass of data, and not the original 
observations themselves, is more nearly the fundamental data "unit” 
within our solution. It was shown that the scaling factor represents 
f= 1/n, where n=50, and approximates reasonably well the average number 
of observations found within a pass of electronic and laser tracking 
data in GEM-T1 . For optical tracking, an n=10 was used through a 
rescaling of data weights by a factor of 5; this was applied to properly 
reflect the lower data density found within this observation type. 

Appendix A2 has developed a general mathematical form of the error 
covariance for the GEM-T1 solution parameters which includes effects of 
unmodeled biases within the data. It was shown that the simplified 

_ I 

inverse of the weighted least - squares normal equations (N“ ) can 
approximate the error covariance by downweighting the data to account 
for the unmodeled bias effects. The weighting should account for the 
effect of a bias per pass of tracking data to be most complete. 
However, it is reasonable to statistically account for the biases by 
applying an average weighting factor to each different satellite data 
type. This latter process was achieved somewhat in GEM-T1 to 
approximate the error covariance matrix. It was also shown that the use 
of least squares collocation ( a priori constraint on the size of the 
coefficients using Kaula's rule) which introduces its own bias-like 
influence on the coefficient solution (favoring zero power) is well 
accommodated by the simplified inverse. The process was also shown to 
perform well because the biases behave largely in a random way across 
the solution. 
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Appendix A3 developed an optimized technique for estimating a 
weighting factor for the different data types in the weighted least- 

squares solution. The algorithm indicates that by iterating the data- 

_ 1 

weights to a converged solution, the error covariance matrix (N ) will 
automatically be calibrated. A solution which properly combines the 
subsets of data at their appropriate weights will also be achieved. 
Again, GEM-T1 was shown to approximate this process in the estimation of 
its data weighting factors. 

Some of the benefits and recommendations which can be made on the 
basis of the three parts of Appendix A are: 

1) It was shown that the more recent data has a diminished 
biased behaviour which leads us to believe that older data 
sets can be replaced as new data becomes available; the 
solution will likewise improve. 

2) We are attempting to understand the modeling errors and 
system causes for these biases so that they can be either 
eliminated or specifically recovered within the solution, 
thusly reducing the need for any downweighting of the least- 
squares normal equations. 

3) We will use the methods developed herein as an algorithm for 
optimal data weighting. Future GEM models will be easier to 
develop and less difficult to optimize given the analysis 
performed here. 
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